#192 4th order Runge-Kutta method

[code language=”python”] # 4th order Runge-Kutta
import os
import numpy as np
import matplotlib.pyplot as plt
import math
os.chdir(r’D:\projects\wordpress\ex192′)
os.getcwd()
# initialization
a=0.0
b=10.0
n=10000
ydumb=np.zeros((2),float)
y=np.zeros((2),float)
fReturn=np.zeros((2),float)
k1=np.zeros((2),float)
k2=np.zeros((2),float)
k3=np.zeros((2),float)
k4=np.zeros((2),float)
y[0]=3.0
y[1]=-5.0
t=a
h=(b-a)/n

def f(t,y):
# force function
fReturn[0] = y[1] fReturn[1] = -100.0 *y[0] -2.0*y[1] + 10.0*math.sin(3*t)
return fReturn

def rk4(t,h,n):
k1=[0]*(n)
k2=[0]*(n)
k3=[0]*(n)
k4=[0]*(n)
fR=[0]*(n)
ydumb=[0]*(n)
# return RHS
fR=f(t,y)
for i in range(0,n):
k1[i]=h*fR[i] for i in range(0,n):
ydumb[i]=y[i] + k1[i]/2
k2=h*f(t+h/2,ydumb)
for i in range(0,n):
ydumb[i]=y[i] + k2[i]/2
k3=h*f(t+h/2,ydumb)
for i in range(0,n):
ydumb[i]=y[i] + k3[i]/2
k4=h*f(t+h,ydumb)
for i in range(0,2):
y[i]=y[i] + (k1[i] + 2.0* (k2[i] + k3[i]) + k4[i])/6.0
return y

tarray=[] y1array=[] y2array=[] while (t < b):
if ((t+h) > b):
h = b-t
y=rk4(t,h,2)
t = t+h
tarray.append(t)
y1array.append(y[0])
y2array.append(y[1])

plt.plot(tarray,y1array,color=’r’,label="y1")
plt.plot(tarray,y2array,color=’g’,label="y2")
plt.xlim(a, b)
plt.legend(["4nd Order RK"], loc=1)
plt.xlabel(‘t’, fontsize=17)
plt.ylabel(‘y(t)’, fontsize=17)
plt.legend()
plt.tight_layout()
plt.savefig("example192_1.png", dpi=300)
plt.show()
plt.close()
plt.plot(y1array,y2array)
plt.legend(["4nd Order RK"], loc=1)
plt.xlabel(‘y1’, fontsize=17)
plt.ylabel(‘y2’, fontsize=17)
plt.tight_layout()
plt.savefig("example192_2.png", dpi=300)
plt.show()
plt.close()
[/code]