def f1(x, y1, y2,y3):
df = -0.013y1-1000y1*y2
return df
def f2(x, y1, y2,y3):
df = -2500y2y3
return df
def f3(x,y1,y2,y3):
df=-0.013y1-1000y1y2-2500y2*y3
return df
def RK4(x, y1, y2,y3, h):
“”"
:param x: Initial value of X
:param y1: Initial value of y1
:param y2: Initial value of y2
:param y3: Initial value of y3
:param h: time step
:return: New iterative solution
“”"
xarray, y1array, y2array,y3array = [], [], [], []
while x <= 0.02:
xarray.append(x)
y1array.append(y1)
y2array.append(y2)
y3array.append(y3)
x += h
K_1 = f1(x, y1, y2,y3)
L_1 = f2(x, y1, y2,y3)
M_1 = f3(x, y1, y2,y3)
K_2 = f1(x + h / 2, y1 + h / 2 * K_1, y2 + h / 2 * L_1 , y3 + h/2 * M_1)