第四十六篇 常微分方程的一步法

之所以称为一步方法，是因为只需要有关前一步的信息就可以在下一步生成解。这使得一步法在计算机程序中相对容易实现。正如许多数值方法的典型，每一步做的工作越多，通常获得的精度就越大。需要在增加每个步骤的工作和减少范围的步骤数之间进行权衡。

一步法主要分为下面几种类型

欧拉法

修正欧拉法

中点法

4阶龙格库塔法

详情可见上篇文章常微分方程编程基础

程序如下

#常微分方程的一步法import numpy as npitype=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 f71print('常微分方程的一步法')if itype==1:print('**欧拉法**')print('xy(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+helif 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+helif 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+helif 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

终端输出结果如下：