Runge Kutta 4th order Python

Question

I am trying to solve this equation using Runge Kutta 4th order:

applying d2Q/dt2=F(y,x,v) and dQ/dt=u Q=y in my program.

I try to run the code but i get this error:

Traceback (most recent call last):
  File "C:\Users\Egw\Desktop\Analysh\Askhsh1\asdasda.py", line 28, in <module>
    k1 = F(y, u, x)  #(x, v, t)
  File "C:\Users\Egw\Desktop\Analysh\Askhsh1\asdasda.py", line 13, in F
    return ((Vo/L -(R0/L)*u -(R1/L)*u**3 - y*(1/L*C)))
OverflowError: (34, 'Result too large')

I tried using the decimal library but I still couldnt make it work properly.I might have not used it properly tho.

My code is this one:

import numpy as np
from math import pi
from numpy import arange
from matplotlib.pyplot import plot, show
#parameters
R0 = 200
R1 = 250
L = 15
h = 0.002
Vo=1000
C=4.2*10**(-6)
t=0.93

def F(y, u, x):
    return ((Vo/L -(R0/L)*u -(R1/L)*u**3 - y*(1/L*C)))


xpoints = arange(0,t,h)
ypoints = []
upoints = []

y = 0.0
u = Vo/L

for x in xpoints:
    ypoints.append(y)
    upoints.append(u)

    m1 = u
    k1 = F(y, u, x)  #(x, v, t)

    m2 = h*(u + 0.5*k1)
    k2 = (h*F(y+0.5*m1, u+0.5*k1, x+0.5*h))

    m3 = h*(u + 0.5*k2)
    k3 = h*F(y+0.5*m2, u+0.5*k2, x+0.5*h)

    m4 = h*(u + k3)
    k4 = h*F(y+m3, u+k3, x+h)

    y += (m1 + 2*m2 + 2*m3 + m4)/6
    u += (k1 + 2*k2 + 2*k3 + k4)/6

plot(xpoints, upoints)
show()

plot(xpoints, ypoints)
show()

I expected to get the plots of u and y against t.

Answer 1

If I may draw your attention to these 4 lines

    m1 = u
    k1 = F(y, u, x)  #(x, v, t)

    m2 = h*(u + 0.5*k1)
    k2 = (h*F(y+0.5*m1, u+0.5*k1, x+0.5*h))

You should note a fundamental structural difference between the first two lines and the second pair of lines.

You need to multiply with the step size h also in the first pair.

---

The next problem is the step size and the cubic term. It contributes a term of size 3*(R1/L)*u^2 ~ 50*u^2 to the Lipschitz constant. In the original IVP per the question with u=Vo/L ~ 70 this term is of size 2.5e+5. To compensate only that term to stay in the stability region of the method, the step size has to be smaller 1e-5.

In the corrected initial conditions with u=0 at the start the velocity u remains below 0.001 so the cubic term does not determine stability, this is now governed by the last term contributing a Lipschitz term of 1/sqrt(L*C) ~ 125. The step size for stability is now 0.02, with 0.002 one can expect quantitatively useful results.

Answer 2

Turns out I messed up with the equations I was using for Runge Kutta

The correct code is the following:

import numpy as np
from math import pi
from numpy import arange
from matplotlib.pyplot import plot, show
#parameters
R0 = 200
R1 = 250
L = 15
h = 0.002
Vo=1000
C=4.2*10**(-6)
t0=0
#dz/dz
def G(x,y,z):
    return Vo/L -(R0/L)*z -(R1/L)*z**3 - y/(L*C)
#dy/dx
def F(x,y,z):
        return z



t = np.arange(t0, 0.93, h)
x = np.zeros(len(t))
y = np.zeros(len(t))
z = np.zeros(len(t))

y[0] = 0.0
z[0] = 0

for i in range(1, len(t)):

        k0=h*F(x[i-1],y[i-1],z[i-1])
        l0=h*G(x[i-1],y[i-1],z[i-1])
        k1=h*F(x[i-1]+h*0.5,y[i-1]+k0*0.5,z[i-1]+l0*0.5)
        l1=h*G(x[i-1]+h*0.5,y[i-1]+k0*0.5,z[i-1]+l0*0.5)
        k2=h*F(x[i-1]+h*0.5,y[i-1]+k1*0.5,z[i-1]+l1*0.5)
        l2=h*G(x[i-1]+h*0.5,y[i-1]+k1*0.5,z[i-1]+l1*0.5)
        k3=h*F(x[i-1]+h,y[i-1]+k2,z[i-1]+l2)
        l3 = h * G(x[i - 1] + h, y[i - 1] + k2, z[i - 1] + l2)

        y[i]=y[i-1]+(k0+2*k1+2*k2+k3)/6
        z[i] = z[i - 1] + (l0 + 2 * l1 + 2 * l2 + l3) / 6
Q=y 
I=z 
plot(t, Q)
show()

plot(t, I)
show()

Answer 3

You can use decimal libary for more precision (handle more digits), but it's kind of annoying every value should be the same class (decimal.Decimal).

For example:

import numpy as np
from math import pi
from numpy import arange
from matplotlib.pyplot import plot, show

# Import decimal.Decimal as D
import decimal
from decimal import Decimal as D

# Precision
decimal.getcontext().prec = 10_000_000

#parameters

# Every value should be D class (decimal.Decimal class)
R0 = D(200)
R1 = D(250)
L = D(15)
h = D(0.002)
Vo = D(1000)
C = D(4.2*10**(-6))
t = D(0.93)

def F(y, u, x):
    # Decomposed for use D
    a = D(Vo/L)
    b = D(-(R0/L)*u)
    c = D(-(R1/L)*u**D(3))
    d = D(-y*(D(1)/L*C))
    return ((a + b + c + d ))


xpoints = arange(0,t,h)
ypoints = []
upoints = []

y = D(0.0)
u = D(Vo/L)

for x in xpoints:
    ypoints.append(y)
    upoints.append(u)

    m1 = u
    k1 = F(y, u, x)  #(x, v, t)

    m2 = (h*(u + D(0.5)*k1))
    k2 = (h*F(y+D(0.5)*m1, u+D(0.5)*k1, x+D(0.5)*h))

    m3 = h*(u + D(0.5)*k2)
    k3 = h*F(y+D(0.5)*m2, u+D(0.5)*k2, x+D(0.5)*h)

    m4 = h*(u + k3)
    k4 = h*F(y+m3, u+k3, x+h)

    y += (m1 + D(2)*m2 + D(2)*m3 + m4)/D(6)
    u += (k1 + D(2)*k2 + D(2)*k3 + k4)/D(6)

plot(xpoints, upoints)
show()

plot(xpoints, ypoints)
show()

But even with ten million of precision I still get an overflow error. Check the components of the formula, their values are way too high. You can increase precision for handle them, but you'll notice it takes time to calculate them.

Answer 4

Problem implementation using scipy.integrate.odeint and scipy.integrate.solve_ivp.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp

# Input data initial conditions
ti = 0.0
tf = 0.5
N  = 100000
h  = (tf-ti)/N

# Initial conditions
u0 = 0.0
Q0 = 0.0

t_span = np.linspace(ti,tf,N)
r0     = np.array([Q0,u0])

# Parameters
R0 = 200
R1 = 250
L  = 15
C  = 4.2*10**(-6)
V0 = 1000

# Systems of First Order Equations

# This function is used with odeint, as specified in the documentation for scipy.integrate.odeint
def f(r,t,R0,R1,L,C,V0):
    Q,u = r
    ode1 = u
    ode2 = -((R0/L)*u)-((R1/L)*u**3)-((1/(L*C))*Q)+(V0/L)
    return np.array([ode1,ode2])

# This function is used in our 4Order Runge-Kutta implementation and in scipy.integrate.solve_ivp
def F(t,r,R0,R1,L,C,V0):
    Q,u = r
    ode1 = u
    ode2 = -((R0/L)*u)-((R1/L)*u**3)-((1/(L*C))*Q)+(V0/L)
    return np.array([ode1,ode2])

# Resolution with oedint
sol_1 =    odeint(f,r0,t_span,args=(R0,R1,L,C,V0))
sol_2 = solve_ivp(fun=F,t_span=(ti,tf), y0=r0, method='LSODA',args=(R0,R1,L,C,V0))

Q_odeint, u_odeint       = sol_1[:,0], sol_1[:,1]
Q_solve_ivp, u_solve_ivp = sol_2.y[0,:], sol_2.y[1,:]

# Figures
plt.figure(figsize=[30.0,10.0])

plt.subplot(3,1,1)
plt.grid(color = 'red',linestyle='--',linewidth=0.4)
plt.plot(t_span,Q_odeint,'r',t_span,u_odeint,'b')
plt.xlabel('t(s)')
plt.ylabel('Q(t), u(t)')

plt.subplot(3,1,2)
plt.plot(sol_2.t,Q_solve_ivp,'g',sol_2.t,u_solve_ivp,'y')
plt.grid(color = 'yellow',linestyle='--',linewidth=0.4)
plt.xlabel('t(s)')
plt.ylabel('Q(t), u(t)')

plt.subplot(3,1,3)
plt.plot(Q_solve_ivp,u_solve_ivp,'green')
plt.grid(color = 'yellow',linestyle='--',linewidth=0.4)
plt.xlabel('Q(t)')
plt.ylabel('u(t)')

plt.show()

Runge-Kutta 4th

# Code development of Runge-Kutta 4 Order
# Parameters
R0 = 200
R1 = 250
L  = 15
C  = 4.2*10**(-6)
V0 = 1000

# Input data initial conditions #
ti = 0.0
tf = 0.5
N  = 100000
h  = (tf-ti)/N

# Initial conditions
u0 = 0.0
Q0 = 0.0

# First order ordinary differential equations
def f1(t,Q,u):
    return u

def f2(t,Q,u):
    return -((R0/L)*u)-((R1/L)*u**3)-((1/(L*C))*Q)+(V0/L)

t = np.zeros(N); Q = np.zeros(N); u = np.zeros(N)
t[0] = ti
Q[0] = Q0
u[0] = u0

for i in range(0,N-1,1):

    k1 = h*f1(t[i],Q[i],u[i])
    l1 = h*f2(t[i],Q[i],u[i])

    k2 = h*f1(t[i]+(h/2),Q[i]+(k1/2),u[i]+(l1/2))
    l2 = h*f2(t[i]+(h/2),Q[i]+(k1/2),u[i]+(l1/2))

    k3 = h*f1(t[i]+(h/2),Q[i]+(k2/2),u[i]+(l2/2))
    l3 = h*f2(t[i]+(h/2),Q[i]+(k2/2),u[i]+(l2/2))

    k4 = h*f1(t[i]+h,Q[i]+k3,u[i]+l3)
    l4 = h*f2(t[i]+h,Q[i]+k3,u[i]+l3)

    Q[i+1] = Q[i] + ((k1+2*k2+2*k3+k4)/6)
    u[i+1] = u[i] + ((l1+2*l2+2*l3+l4)/6)
    t[i+1] = t[i] + h

plt.figure(figsize=[20.0,10.0])
plt.subplot(1,2,1)
plt.plot(t,Q_solve_ivp,'r',t,Q_odeint,'y',t,Q,'b')
plt.grid(color = 'yellow',linestyle='--',linewidth=0.4)
plt.xlabel('t(s)')
plt.ylabel(r'$Q(t)_{Odeint}$, $Q(t)_{RK4}$')

plt.subplot(1,2,2)
plt.plot(t,Q_solve_ivp,'g',t,Q_odeint,'y',t,Q,'b')
plt.grid(color = 'yellow',linestyle='--',linewidth=0.4)
plt.xlabel('t(s)')
plt.ylabel(r'$Q(t)_{solve_ivp}$, $Q(t)_{RK4}$')