Runge Kutta Fehlberg goes crazy when using adaptive h in python

Question

So basically what i'm trying to do here is try to define a function that integrates for runge kutta. The most "physical" aspects of it i have already controlled for (it's not a problem of the values of the functions , or the physical parameters (initial conditions, etc), since i have controlled for those and they don't change the result) . The problem i'm solvin is an ODE of second order (harmonic oscilator),

x''-μ(1-x^2 )x'+x=0 + CI x(0)=2, x'(0)=0, for t in [0,25], and μ=4

which i have separated into two first-order ODES;

μ*(1-x**2)*v-x = v'
v = x'

and we have inital conditions x0,v0,t0 and a interval on which to integrate t0,tf, and a step size h with which to integrate. ε is my error tolerance. The error I have is that when i activate the adaptivetoggle, so that i should use an adaptive h, x just goes to a value and v goes to 0, instead of oscilating like a harmonic oscilator should. I suspect the problem is only about that bit of code, because when i deactivate the adaptive toggle , everything runs just fine. I am not sure of what is happening, my values should be oscilating, and small, and instead the error just goes to 0 ( it shouldn't, i think ) , v does too and x tends to something instead of oscilating.

The code i'm running is:

def rk452D(v, f, x0, v0, t0, tf, ε, h, adaptivetoggle):
x = x0; t = t0 
valuesx = []
valuesv = []
while t<tf :
#We define the runge kutta functions on which to iterate while t is in [t0,tf], and they are of the kind (t,x,v)

    f1v = f(t,               x,                                                                                    v0                                                                                       )
    f1  = v(t,               x,                                                                                    v0                                                                                       )
f2v = f(t + (1/4)*h,     x + (1/4)*f1,                                                                         v0 + (1/4)*f1v                                                                           )
f2  = v(t + (1/4)*h,     x + (1/4)*f1,                                                                         v0 + (1/4)*f1v                                                                           )

f3v = f(t + (3/8)*h,     x + (3/32)*f1      + (9/32)*f2  ,                                                     v0 + (3/32)*f1v      + (9/32)*f2v                                                        )
f3  = v(t + (3/8)*h,     x + (3/32)*f1      + (9/32)*f2  ,                                                     v0 + (3/32)*f1v      + (9/32)*f2v                                                        )

f4v = f(t + (12/13)*h,   x + (1932/2197)*f1 - (7200/2197)*f2 + (7296/2197)*f3,                                 v0 + (1932/2197)*f1v - (7200/2197)*f2v + (7296/2197)*f3v                                 )
f4  = v(t + (12/13)*h,   x + (1932/2197)*f1 - (7200/2197)*f2 + (7296/2197)*f3,                                 v0 + (1932/2197)*f1v - (7200/2197)*f2v + (7296/2197)*f3v                                 )

f5v = f(t + h,           x + (439/216)*f1   - 8*f2           + (3680/513)*f3  - (845/4104)*f4,                 v0 + (439/216)*f1v   - 8*f2v           + (3680/513)*f3v  - (845/4104)*f4v                )
f5  = v(t + h,           x + (439/216)*f1   - 8*f2           + (3680/513)*f3  - (845/4104)*f4,                 v0 + (439/216)*f1v   - 8*f2v           + (3680/513)*f3v  - (845/4104)*f4v                )

f6v = f(t + h/2,         x - (8/27)*f1      + 2*f2           - (3544/2565)*f3 + (1859/4104)*f4 - (11/40)*f5,   v0 - (8/27)*f1v      + 2*f2v           - (3544/2565)*f3v + (1859/4104)*f4v - (11/40)*f5v )
f6  = v(t + h/2,         x - (8/27)*f1      + 2*f2           - (3544/2565)*f3 + (1859/4104)*f4 - (11/40)*f5,   v0 - (8/27)*f1v      + 2*f2v           - (3544/2565)*f3v + (1859/4104)*f4v - (11/40)*f5v )


#Now we calculate the positions and velocities, for the fourth order runge kutta aproximation. Commented we have the fifth order approxiation results, which we use to estimate the error ( Error = abs(order5approximation-order4approximation)
    x4 = x + h((25/216)f1  + (1408/2565)f3   + (2197/4104)f4     -(1/5)f5              )
    #x5 = x + h((16/135)f1  + (6656/12825)f3  + (28561/56430)f4   -(9/50)f5  +(2/55)*f6 )
v4 = v0 + h*((25/216)*f1v + (1408/2565)*f3v  + (2197/4104)*f4v    -(1/5)*f5v             )
#v5 = v0 + h*((16/135)*f1v + (6656/12825)*f3v + (28561/56430)*f4v  -(9/50)*f5v +(2/55)*f6v)        


#If we want to use an adaptive h for our calculations,
    if adaptivetoggle == True :
       #We calculate error in x using fs, eror in v using fvs, and we take the smaller hnew of the two 
       Errorx = abs((1/360)f1  - (128/4275)f3  - (2197/75240)f4  + (1/50)f5  + (2/55)f6  )
       Errorv = abs((1/360)f1v - (128/4275)f3v - (2197/75240)f4v + (1/50)f5v + (2/55)f6v )
   hnewx  = 0.9*h*((ε)/Errorx)**(0.25)
   hnewv  = 0.9*h*((ε)/Errorv)**(0.25)

   hnew = min(hnewx,hnewv)

   print(hnew)
   #if hnewx &lt; hnewv:
    #   hnew = hnewx
   #if hnewx &gt; hnewv:
     #  hnew = hnewv

   #After calculating hnew, we compare it to the integration step h we have used, to decide if we can keep our calculation or we have to repeat it using hnew.

   if hnew &gt;= h :
      #we increment the loop,and take hnew for the next loop
      t += h
      h = hnew
      x = x4
      v0 = v4
      valuesx.append(x4)
      valuesv.append(v4)

   elif hnew &lt; h: 
      h  = hnew #we don't increment t , the loop variable, when we repeat the integration step using the new h integration step

else :#if we don't want an adaptive h ( this works just fine)
   valuesx.append(x4)
   valuesv.append(v4)
   x = x4
   v0 = v4
   t+=h #increment the loop

return valuesx, valuesv
#Then we implement the function
#We define the two functions ( of the ODEs) to feed the RK solver
def f(t,x,v):
    return μ(1-x2)v-x
def v(t,x,v):
    return v
#we feed it the parameters
μ  = 4; x0 = 2;  v0 = 0; t0 = 0; tf = 25; h  = 1; ε  = 10**-3 
adaptivetoggle = True
solution = rk452D(v, f, x0, v0, t0, tf, ε, h, adaptivetoggle)
print(solution)

Lutz Lehmann · Answer 1 · 2020-12-08T09:45:58.653

The error is so glaringly blatant, so blending and ubiquitous that it has become almost invisible. Straight-up not seeing the forest for the trees.

In the computation of the slopes f1, f1v, ... you do not use the step size h. You can see the effect of that if you enlarge the debug print command to also print t,h and Errorx, Errorv.

For better readability, single-point-of-failure and slightly better efficiency I'd suggest changing to the format

    ttmp = t + (12/13)*h
    xtmp = x  + h*((1932/2197)*f1  - (7200/2197)*f2  + (7296/2197)*f3 )
    vtmp = v0 + h*((1932/2197)*f1v - (7200/2197)*f2v + (7296/2197)*f3v)
    f4v, f4 = f(ttmp, xtmp, vtmp ), v(ttmp, xtmp, vtmp )

(or use a general vectorized version, without treating the components separately).

Then the initial step size h=1 is still too large as the controller regulates down to h=0.05, periodically raising up to h=0.6.

For faster testing, the modified slope computation is

        f1v = f(t,  x,  v0 )
        f1  = v(t,  x,  v0 )
    f2v = f(t + (1/4)*h,  x + h*(1/4)*f1,  v0 + h*(1/4)*f1v )
    f2  = v(t + (1/4)*h,  x + h*(1/4)*f1,  v0 + h*(1/4)*f1v )

    tstg = t + (3/8)*h   # stg like stage
    xstg = x  + h*((3/32)*f1  + (9/32)*f2 )
    vstg = v0 + h*((3/32)*f1v + (9/32)*f2v)
    f3v, f3 = f(tstg, xstg,  vstg), v(tstg, xstg,  vstg)

    tstg = t + (12/13)*h
    xstg = x  + h*((1932/2197)*f1  - (7200/2197)*f2  + (7296/2197)*f3 )
    vstg = v0 + h*((1932/2197)*f1v - (7200/2197)*f2v + (7296/2197)*f3v)
    f4v, f4 = f(tstg, xstg,  vstg), v(tstg, xstg,  vstg)

    tstg = t + h
    xstg = x  + h*((439/216)*f1  - 8*f2  + (3680/513)*f3  - (845/4104)*f4 )
    vstg = v0 + h*((439/216)*f1v - 8*f2v + (3680/513)*f3v - (845/4104)*f4v)
    f5v, f5 = f(tstg, xstg,  vstg), v(tstg, xstg,  vstg)

    tstg = t + h/2
    xstg = x  + h*(-(8/27)*f1  + 2*f2  - (3544/2565)*f3  + (1859/4104)*f4  - (11/40)*f5 )
    vstg = v0 + h*(-(8/27)*f1v + 2*f2v - (3544/2565)*f3v + (1859/4104)*f4v - (11/40)*f5v)
    f6v, f6 = f(tstg, xstg,  vstg), v(tstg, xstg,  vstg)

One might change the acceptance condition to hnew >= 0.95*h or hnew >= 0.85*h to avoid too many step repetitions, with the factor 0.9 in the computation of hnew there is enough safety to allow that.

May thanks!!! What do i mean i don't use the step size h? you mean it should be fi=h*f()? Isn't it good to use just an h at the calculation of x,v, since every function has a h multiplying, and so i can just put it into x,v? — L.Change, Dec 08 '20 at 10:01
I am going to test it, i dont fully understand your second code tho ( pardon me, i'm not too good) , you sseem to calculate x, v for rk2,rk3,rk4? or why do you calculate x,v every time you calculate a function, just to test ? — L.Change, Dec 08 '20 at 10:03
Compare my computation of f2 with your original code. Yes, in each computation of the x value the update has to be scaled by h, either in front of the function or inside. In x4 you used the "inside" variant, that is, the f values are not scaled, so I extended this to all occurrences. // I just took the linear combinations out of the function arguments, using temporary variables to hold their values. The computation stays the same (apart from the factor h). — Lutz Lehmann, Dec 08 '20 at 10:09
Yeah i ran the code with the f(h*f1+etc etc) as you said, and it is indeed running at least logically, it oscilates. I have to plot this and everything, but dude, you are awesome and thank you very much. I spent 7+hours putting prints everywhere to test where it went wrong, and also controlling for any possibilities (maybe the parameters, etc etc). Thank you so very much man. I'm going to read a bit better your code, to try and understand what you're doing there. Thanks! — L.Change, Dec 08 '20 at 10:13
You should always test your self-programmed methods for basic assumptions, for instance like in stackoverflow: DoPri45 with fixed step or a 3rd order method and the chaotic Lorenz attractor — Lutz Lehmann, Dec 08 '20 at 14:40
yeah you were totally right....now i'm implementing it in 1D, and have corrected that. I programmed some plots, and corrected some bugs and errors i dod on those, but I have a very weird thing were the adaptive step size method is giving me a worse approximation, and also the program says the error is always zero in the plots, which it shouldn't be-it should oscilate, and after some hours i'm still unable to see the mistake, I posted a question about it in the meantime to see if anyone smarter than me can figure it out faster... — L.Change, Dec 08 '20 at 20:26

Mason · Answer 2 · 2020-12-08T18:10:30.573

1

The code doesn't run as it is now. I think your formula for the error is wrong. Judging by the $4$th root you took, the tolerance $\varepsilon$ is meant to be for error per time since error per time is $O(h^4)$ for RK4. So you should use Error/h instead of Error in your formula for hnew. Also your criteria for when to accept the approximation is harsh (it may not accept an approximation even if it was within the tolerance). You could just check whether Error <= $\varepsilon$.

Also you could take the maximum of the errors before computing the hnew rather than computing two hnew and taking a minimum. This error would be the $l_{\infty}$ norm of (RK4 - RK5)/h, and the hnew formula is based on the fact that the error is $O(h^4)$.

edited Dec 08 '20 at 18:10

answered Dec 08 '20 at 05:23

Mason

10,415

As in the computation of the error, the multiplication by $h$ is omitted, there is no need for a division by $h$ afterwards. – Lutz Lehmann Dec 08 '20 at 05:53
Thank you! The code should run though, but give very bad results. In my error calculation i take the absolute value of the order 5 - order 4 solutions , and since they are h*something, i just ommited the h and wrote the coefficients directly instead of leaving it and then using h in my formula for hnew. – L.Change Dec 08 '20 at 09:09
I'm going to try and chech error/h, and also change the election from the h to the error, although, the problem is that the error itself , when i only chang the adaptivetoggle to false, explodes – L.Change Dec 08 '20 at 09:10

Runge Kutta Fehlberg goes crazy when using adaptive h in python

2 Answers2