One of my first blogs was applying Parker-Sochacki to circular motion where dealt with the following couple Ordinary Differential Equations:

x'(t) = y(t) ; x(0) = 1

y'(t) = -x(t) ; y(0) = 0

In that blog, I compared solving this with Forward Euler, Runge-Kutta Order 2, Runge-Kutta Order 3, Runge-Kutta 4, and different orders of Parker-Sochacki. I showed implementation in Matlab and Houdini. Inside Houdini, I used VOPs.

I wanted to show this using C++ like code in the Attribute Wrangle node of Houdini.

Recently, I started using the Attribute Wrangle node in Houdini. The Attribute Wrangle node uses VEX code, which is similar to what you see in C++. Most of the code I have shown is in Matlab, which can be confusing because the formulas I show are different than Matlab. Matlab starts at index 1 instead of 0.

I am reshowing the circle implementation here: how to compute the coefficients, how to generate the Maclaurin Series, and a screenshot of the VEX code for you to reference if you are using Parker-Sochacki in a language like C++ or Python.

x_i+1 = -y_i / (i+1)

y_i+1 = x_i / (i+1)

x(t) = x_0t^0 + x_1t^1 + x_2t^2 + x_3t^3 + … + x_nt^n

y(t) = y_0t^0 + y_1t^1 + y_2t^2 + y_3t^3 + … + y_nt^n

n is a positive integer of your choosing. Generally speaking, a higher n value leads to a higher accuracy. This does increase computation time. t is the time step, which is also something you can pick. I chose 1/24 because this is Houdini’s default frame rate. Decreasing the time step increases the accuracy but the whole point of Parker-Sochacki(with what I am doing anyway) is increasing the order so I can take a larger time step saving time overall.

Here is the code. I am not claiming this is the optimal way to do this. This is a way to do this. The arrays were prebuilt.

In case you can’t see it:

————————————————————-

int order = 6;

float x = @P.x;

float y = @P.y;

float xcoeff[] = {0,0,0,0,0};

float ycoeff[] = {0,0,0,0,0};

xcoeff[0] = @P.x;

ycoeff[0] = @P.y;

float special_time[] = {1,0.041666667, 0.0017361111, 0.00007233796296296, 0.000003014081790123456, 0.000000125586741255144032921810699588477366};

for (int i = 0; i < order-1; i++)

{

xcoeff[i+1] = -ycoeff[i] / (i+1);

ycoeff[i+1] = xcoeff[i] / (i+1);

}

float power_x = 0;

float power_y = 0;

for (int i = 0; i < order-1; i++)

{

power_x += xcoeff[i]*special_time[i];

power_y += ycoeff[i]*special_time[i];

}

@P.x = power_x;

@P.y = power_y;

————————————————————-

@P.x and @P.y are already defined in Houdini to update the x position and y position.