1

I'm doing a homework assignment in which I have to model the motion of an aircraft with a 4th order RK approximation. However, for some reason I keep getting results that indicate that all of my K values are zero when I iterate.

my main loop is as follows:

for(i=0;i<duration;i++){
    t[i]=(i)*delta_t;
    rho[i]=1.225*(exp((-0.0001)*z[i]));

    Kx1=f(i, V_x[i], V_z[i], rho[i], z[i], m, S, Cl);
    Kz1=g(i, V_x[i], V_z[i], rho[i], z[i], m, S, Cl, G);
    Kx2=f(i, (V_x[i]+((Kx1*h)/2)), (V_z[i]+((Kz1*h)/2)), rho[i], z[i], m, S, Cl);
    Kz2=g(i, (V_x[i]+((Kx1*h)/2)), (V_z[i]+((Kz1*h)/2)), rho[i], z[i], m, S, Cl, G);
    Kx3=f(i, (V_x[i]+((Kx2*h)/2)), (V_z[i]+((Kz2*h)/2)), rho[i], z[i], m, S, Cl);
    Kz3=g(i, (V_x[i]+((Kx2*h)/2)), (V_z[i]+((Kz2*h)/2)), rho[i], z[i], m, S, Cl, G);
    Kx4=f(i, (V_x[i]+((Kx3*h)/2)), (V_z[i]+((Kz3*h)/2)), rho[i], z[i], m, S, Cl);
    Kz4=g(i, (V_x[i]+((Kx3*h)/2)), (V_z[i]+((Kz3*h)/2)), rho[i], z[i], m, S, Cl, G);

    V_x[i+1]=V_x[i]+(delta_t/6)*(Kx1+2*Kx2+2*Kx3+Kx4);
    x[i+1]=x[i]+(V_x[i+1])*delta_t;
    V_z[i+1]=V_z[i]+(delta_t/6)*(Kz1+2*Kz2+2*Kz3+Kz4);
    z[i+1]=z[i]+(V_z[i+1])*delta_t;

    t[0]=0;
    x[0]=0;
    z[0]=4800;
    V_z[0]=0;
    V_x[0]=210;
}

with V_z, V_x, x, z, and rho all being 1xduration arrays (I have duration defined as 400). f and g are two functions that return acceleration in x and z directions respectively. I've tested both externally and by hand and they work. I reinforced the boundary conditions in the loop (I've written them outside as well) because I was getting weird results otherwise, but that's beside the point of what I'm having trouble with: my K values all seem to be ending up equal to zero because V_x stays at 210 for each iteration, V_z stays at 0, and obviously then x and z don't change.

I created a quick test code in a separate file to calculate the K values without iterating just to make sure the algorithm is sound (all the same BCs, constants, etc., just no loop):

    Kx1=f(i, V_x[0], V_z[0], rho[0], z[0], m, S, Cl);
    Kz1=g(i, V_x[0], V_z[0], rho[0], z[0], m, S, Cl, G);
    Kx2=f(i, (V_x[0]+((Kx1*h)/2)), (V_z[0]+((Kz1*h)/2)), rho[0], z[0], m, S, Cl);
    Kz2=g(i, (V_x[0]+((Kx1*h)/2)), (V_z[0]+((Kz1*h)/2)), rho[0], z[0], m, S, Cl, G);
    Kx3=f(i, (V_x[0]+((Kx2*h)/2)), (V_z[0]+((Kz2*h)/2)), rho[0], z[0], m, S, Cl);
    Kz3=g(i, (V_x[0]+((Kx2*h)/2)), (V_z[0]+((Kz2*h)/2)), rho[0], z[0], m, S, Cl, G);
    Kx4=f(i, (V_x[0]+((Kx3*h)/2)), (V_z[0]+((Kz3*h)/2)), rho[0], z[0], m, S, Cl);
    Kz4=g(i, (V_x[0]+((Kx3*h)/2)), (V_z[0]+((Kz3*h)/2)), rho[0], z[0], m, S, Cl, G);

This returns the proper K values, so I'm at a complete loss as to why my V_x[i] and V_x[i] aren't changing in the loop.

Jacob Shapiro
  • 11
  • 1
  • 1
    never mind lol, I had delta_t as an integer so multiplying by it yielded an integer and my k values rounded to 0 when calculating V_x[i+1] and V_z[i+1]. Changed it to a double and the program works perfectly. – Jacob Shapiro Apr 18 '21 at 18:59
  • Good that you found the result. May I ask why you are reinitializing the variables within the for loop? Just curious. In fact, if you are not using the array for further processing (like say restart from some other timestep,) you are also losing speed by using an array. – sudhish Apr 18 '21 at 19:12
  • I was having issues with an earlier code that I solved by re-stating the boundary conditions in the loop and now I do it reflexively since it inexplicably solved an issue with an earlier version of this program. I've removed them from the loop now and it still works exactly the same. As for why I am using an array, the assignment calls for it. Presumably if I were doing this for a professional application I would need to store those values to use them later. – Jacob Shapiro Apr 18 '21 at 20:10
  • You are using `h` and `delta_t` for the same role, the step size. Do they really have always the same value? You are reducing the order of the method to possibly 1 by not including the `z` integration properly. You can test this by computing the differences for an integration with 1/2 and 2 times the step size. The factor between the errors should be 16, if it is less then the order of the method is not 4. – Lutz Lehmann Apr 19 '21 at 10:49

0 Answers0