I created a differential equation solver (Runge-Kutta 4th order method) in Python. Than I decided to check its results by setting the parameter mu to 0 and looking at the numeric solution that was returned by it.
The problem is, I know that this solution should give an stable oscillation as a result, but instead I get a diverging solution.
The code is presented below. I tried solving this problem (rounding errors from floating point precision) by using numpy float128 data type. But the solver keeps giving me the wrong answer.
The code is:
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
def f(t,x,v):
f = -k/m*x-mu/m*v
return(f)
def g(t,x,v):
g = v
return(g)
def srunge4(t,x,v,dt):
k1 = f(t,x,v)
l1 = g(t,x,v)
k2 = f(t+dt/2, x+k1*dt/2, v+l1*dt/2)
l2 = g(t+dt/2, x+k1*dt/2, v+l1*dt/2)
k3 = f(t+dt/2, x+k2*dt/2, v+l2*dt/2)
l3 = g(t+dt/2, x+k2*dt/2, v+l2*dt/2)
k4 = f(t+dt/2, x+k3*dt, v+l3*dt)
l4 = g(t+dt/2, x+k3*dt, v+l3*dt)
v = v + dt/6*(k1+2*k2+2*k3+k4)
x = x + dt/6*(l1+2*l2+2*l3+l4)
t = t + dt
return([t,x,v])
mu = np.float128(0.00); k = np.float128(0.1); m = np.float128(6)
x0 = np.float128(5); v0 = np.float128(-10)
t0 = np.float128(0); tf = np.float128(1000); dt = np.float128(0.05)
def sedol(t, x, v, tf, dt):
sol = np.array([[t,x,v]], dtype='float128')
while sol[-1][0]<=tf:
t,x,v = srunge4(t,x,v,dt)
sol=np.append(sol,np.float128([[t,x,v]]),axis=0)
sol = pd.DataFrame(data=sol, columns=['t','x','v'])
return(sol)
ft_runge = sedol(t0, x0, v0, tf, dt=0.1)
plt.close("all")
graf1 = plt.plot(ft_runge.iloc[:,0],ft_runge.iloc[:,1],'b')
plt.show()
Am I using numpy float128 in a wrong way?
You are mixing in srunge4 the association of k and l to x and v. Per the function association and the final summation, the association should be (v,f,k) and (x,g,l). This has to be obeyed in the updates of the stages of the method.
In stage 4, it should be t+dt in the first argument. However, as t is not used in the derivative computation, this error has no consequence here.
Also, you are destroying the float128 computation if you set one parameter to a float in the default float64 type in dt=0.1.
The code with these corrections and some further simplifications is
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
mu = np.float128(0.00); k = np.float128(0.1); m = np.float128(6)
x0 = np.float128(5); v0 = np.float128(-10)
t0 = np.float128(0); tf = np.float128(1000); dt = np.float128(0.05)
def f(t,x,v): return -(k*x+mu*v)/m
def g(t,x,v): return v
def srunge4(t,x,v,dt): # should be skutta4, Wilhelm Kutta gave this method in 1901
k1, l1 = (fg(t,x,v) for fg in (f,g))
# here is the essential corrections, x->l, v->k
k2, l2 = (fg(t+dt/2, x+l1*dt/2, v+k1*dt/2) for fg in (f,g))
k3, l3 = (fg(t+dt/2, x+l2*dt/2, v+k2*dt/2) for fg in (f,g))
k4, l4 = (fg(t+dt , x+l3*dt , v+k3*dt ) for fg in (f,g))
v = v + dt/6*(k1+2*k2+2*k3+k4)
x = x + dt/6*(l1+2*l2+2*l3+l4)
t = t + dt
return([t,x,v])
def sedol(t, x, v, tf, dt):
sol = [[t,x,v]]
while t<=tf:
t,x,v = srunge4(t,x,v,dt)
sol.append([t,x,v])
sol = pd.DataFrame(data=np.asarray(sol) , columns=['t','x','v'])
return(sol)
ft_runge = sedol(t0, x0, v0, tf, dt=2*dt)
plt.close("all")
fig,ax = plt.subplots(1,3)
ft_runge.plot(x='t', y='x', ax=ax[0])
ft_runge.plot(x='t', y='v', ax=ax[1])
ft_runge.plot.scatter(x='x', y='v', s=1, ax=ax[2])
plt.show()
It produces the expected ellipse without visually recognizable changes in the amplitudes.