数值分析—四阶龙格库塔python复现

python入门，复现《数值分析》中的欧拉、龙格库塔，这里以经典四阶龙格库塔来示例。

四阶龙格库塔公式

y_(n+1)=y_n+h/6 (K_1+2K_2+2K_3+K_4 )
K_1=f(x_n,y_n)
K_2=f(x_n+h/2 ,y_n+〖hK〗_1/2)
K_3=f(x_n+h/2 ,y_n+〖hK〗_2/2）
K_4=f(x_n+h ,y_n+hK_3)

python代码

下面是python代码

#自行更改微分方程
def f(p,q):
return (-1)*q-p*q**2

#把公式拆分更简洁
def zsj(m,n,h):
K1=f(m,n)
K2=f(m+h*0.5,n+h*0.5*K1)
K3=f(m+h*0.5,n+h*0.5*K2)
K4=f(m+h,n+h*K3)
return K1+2*K2+2*K3+K4

#f=微分方程；求解区间(a,b)；ya=y(0),h=步长
def LB(f,a,b,ya,h):
xi=a
while xi<=b:
if xi==0:
yi=ya
else:
yi = yi+h/6*zsj(xi-0.2,yi,h)
print('y(%s)= %s'%(xi,yi))
xi=xi+h

print(LB(f,0,1,1,0.2))

运行结果