第四十六篇 常微分方程的一步法
之所以称为一步方法,是因为只需要有关前一步的信息就可以在下一步生成解。这使得一步法在计算机程序中相对容易实现。正如许多数值方法的典型,每一步做的工作越多,通常获得的精度就越大。需要在增加每个步骤的工作和减少范围的步骤数之间进行权衡。
一步法主要分为下面几种类型
欧拉法
修正欧拉法
中点法
4阶龙格库塔法
详情可见上篇文章常微分方程编程基础
程序如下
#常微分方程的一步法
import numpy as np
itype=4;n=1;nsteps=5;h=0.1
k0=np.zeros((n))
k1=np.zeros((n))
k2=np.zeros((n))
k3=np.zeros((n))
y=np.zeros((n))
x=0.0;y[:]=(1.0)
def f71(x,y):n=y.shape[0]f71=np.zeros((n))f71[0]=(x+y[0])**2return f71
print('常微分方程的一步法')
if itype==1:print('**欧拉法**')print('x y(i) , i=1',n)for j in range(0,nsteps+1):print('{:9.5e}'.format(x),end=' ')for i in range(1,n+1):print('{:9.5e}'.format(y[i-1]))k0=h*f71(x,y);y=y+k0;x=x+h
elif itype==2:print('**修正欧拉法**')print('x y(i) , i=1',n)for j in range(0,nsteps+1):print('{:9.5e}'.format(x),end=' ')for i in range(1,n+1):print('{:9.5e}'.format(y[i-1]))k0=h*f71(x,y);k1=h*f71(x+h,y+k0);y=y+(k0+k1)/2.0;x=x+h
elif itype==3:print('**中点法**')print('x y(i) , i=1',n)for j in range(0,nsteps+1):print('{:9.5e}'.format(x),end=' ')for i in range(1,n+1):print('{:9.5e}'.format(y[i-1]))k0=h*f71(x,y);k1=h*f71(x+h/2.0,y+k0/2.0);y=y+k1;x=x+h
elif itype==4:print('**四阶龙格库塔法**')print('x y(i) , i=1',n)for j in range(0,nsteps+1):print('{:9.5e}'.format(x),end=' ')for i in range(1,n+1):print('{:9.5e}'.format(y[i-1]))k0=h*f71(x,y);k1=h*f71(x+h/2.0,y+k0/2.0)k2=h*f71(x+h/2.0,y+k1/2.0);k3=h*f71(x+h,y+k2)y=y+(k0+2.0*k1+2.0*k2+k3)/6.0;x=x+h
终端输出结果如下: