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 < hnewv:
# hnew = hnewx
#if hnewx > 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 >= 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 < 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-x**2)*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)
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)$.