Search code examples
pythonoderunge-kutta

Runge Kutta 4th order Python

I am trying to solve this equation using Runge Kutta 4th order: enter image description here

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.

Solution

  • 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()