常微分方程初值问题数值解法[完整公式]（Python）

y=np.zeros((n,1))

x=np.zeros((n,1))

y[0] =y0

x[0] = a

for i in range(1,n,1):

x[i]=a+i*h

y[i]= y[i-1]+h*f(x[i-1],y[i-1])

return x,y

def ModEuler(a,b,f,y0,n): #改进Ｅｕｌｅｒ公式

h=np.abs(b-a)/(n-1)

y=np.zeros((n,1))

x=np.zeros((n,1))

y[0] =y0

x[0] = a

for i in range(1,n,1):

x[i]=a+i*h

y[i]= y[i-1]+h*f(x[i-1],y[i-1])

y[i] = y[i-1]+h/2*(f(x[i-1],y[i-1])+f(x[i],y[i]))

return x,y

def Heun(a,b,f,y0,n): #二阶Ｒｕｎｇｅ—Ｋｕｔｔａ方法：Ｈｅｕｎ公式

h=np.abs(b-a)/(n-1)

y=np.zeros((n,1))

x=np.zeros((n,1))

y[0] =y0

x[0] = a

K1,K2=0,0

for i in range(1,n,1):

x[i]=a+i*h

K1 = f(x[i-1],y[i-1])

K2 = f(x[i-1]+2/3h,y[i-1]+2/3h*K1)

y[i] = y[i-1]+h/4*(K1+3*K2)

return x,y

def Ord3Kutta(a,b,f,y0,n): #三阶Kutta公式

h=np.abs(b-a)/(n-1)

y=np.zeros((n,1))

x=np.zeros((n,1))

y[0] =y0

x[0] = a

K1,K2,K3=0,0,0

for i in range(1,n,1):

x[i]=a+i*h

K1 = f(x[i-1],y[i-1])

K2 = f(x[i-1]+1/2h,y[i-1]+1/2h*K1)

K3 = f(x[i-1]+h,y[i-1]-hK1+2h*K2)

y[i] = y[i-1]+h/6*(K1+4*K2+K