1

I am trying to solve a simple 2nd order DE using 4th order Runge-Kutta on C. My code compiles but the results are kinda awkward and I can't find the mistake. I have repeated the calculations multiple times but I always get the same result. The differential equation I am trying to solve is:

enter image description here

and my code is:

#include<stdio.h>
#include<math.h>

int main (){
    int b,t,h,w,m;
    double c1,c2,k1,k2,k3,k4,j1,j2,j3,j4,Y[15],Z[15];
    FILE * fp = fopen("VRL4.csv", "w");
    h=10000.0;
    Y[1]=0.0;
    Z[1]=1.0;

    for(m=1;m<=11;m++){ /*m is an auxiliary variable */
       b=(2.0*m-2.0)*100.0+10;
       c1=0.2222222222*b;   
       c2=0.14814814815*b;
          for(w=1;w<=11;w++){ /* w is an auxiliary variable */

            t=(w-1)*10000+1;

            k1=Z[w];
            j1=-c1*pow(t,-0.33333)*Z[w]-c2*pow(t,-1.33333)*Y[w];

            k2=Z[w]+0.5*j1;
            j2=-c1*pow(t+0.5*h,-0.33333)*(Z[w]+0.5*j1)-c2*pow(t+0.5*h,-1.33333)*(Y[w]+0.5*k1);

            k3=Z[w]+0.5*j2;
            j3=-c1*pow(t+0.5*h,-0.33333)*(Z[w]+0.5*j2)-c2*pow(t+0.5*h,-1.33333)*(Y[w]+0.5*k2);

            k4=Z[w]+j2;
            j4=-c1*pow(t+h,-0.33333)*(Z[w]+j2)-c2*pow(t+h,-1.33333)*(Y[w]+k2);

            Y[w+1]=Y[w]+0.166667*h*(k1+2.0*k2+2*k3+k4);
            Z[w+1]=Z[w]+0.166667*h*(j1+2.0*j2+2*j3+j4);

            fprintf(fp, "%d,%d,%lf\n",b,t,Z[w]);
            fprintf(fp, "%d,%d,%lf\n",b,t,Y[w]);
        }   
    }   

    return 0;
}
Donald Duck
  • 8,409
  • 22
  • 75
  • 99
Suyama87
  • 85
  • 10
  • 2
    Please post the error you are seeing. What are the incorrect results? What Results are you expetcing? etc.. – Garrigan Stafford Aug 02 '18 at 18:45
  • 3
    Using doubles, 0.33333 is not really very close to 1./3, and 0.16667 is not very close to 1./6. – FredK Aug 02 '18 at 19:47
  • 2
    Your step size (h) is 10000. That is massive! Have you tried changing the value to something smaller like 0.1? Also you don't need to store all the Z[w] and Y[w] values you only need the current and next one – Spoonless Aug 03 '18 at 08:23
  • For a 2DE you need to simultaneously solve the both the first and 2nd order derivatives at each step. Check [this example] (https://math.stackexchange.com/questions/721076/help-with-using-the-runge-kutta-4th-order-method-on-a-system-of-2-first-order-od) on maths stack exchange – Spoonless Aug 03 '18 at 08:35
  • @GarriganStafford I don't get an error. My results just jump from 0 (initial condition) to 10^40. I am certainly not expecting such an intense rise. Maybe something more like an oscillation, or a gradual decrease. – Suyama87 Aug 09 '18 at 16:52
  • @Spoonless It's a physics problem and the value of the step size corresponds to years. I can't change it, because then, it wouldn't be the problem I want to solve. – Suyama87 Aug 09 '18 at 16:53
  • I tried changing the initial condition, but I don't get different results. Can you please help me as to what can be so wrong which I happen not to see as well, here? – Suyama87 Aug 09 '18 at 17:00
  • You might be experiencing underflow then – Garrigan Stafford Aug 09 '18 at 17:29
  • @GarriganStafford This is the first time I'm experiencing this problem. Could you please give me some suggestions as to how to handle underflow? – Suyama87 Aug 09 '18 at 17:41
  • @Ine If you want to know the value in 10000 years that is fine, why do you need to estimate that in one step? Regardless of the source of the problem, large steps - especially for high order derivatives - results in instability. – Spoonless Aug 10 '18 at 08:49
  • 1
    @Spoonless Well, first mistake I noticed: When I calculate k4 and j4 I should be using k3 and j3 (and not k2 and j2). But this doesn't really change things. Choosing such a large step was inevitable, since my program crashes otherwise. What bugs me at the moment is that when I plot the solution (Y[w] VS t) the derivative (Z[w]) seems to be negative, while numerically it appears to be positive. – Suyama87 Aug 10 '18 at 15:35
  • The Lipschitz constant is somewhere about `L=1000`. To get about 8 valid digits in the result of RK4, one should one should use `Lh=0.01`, which implies a step size `h=1e-5`, the proposed step size `1e+4` is catastrophically too large. // To ameliorate this use an implicit method with step size control to get the optimal step sizes for method and accuracy target. – Lutz Lehmann Dec 15 '22 at 14:10

0 Answers0