数值分析实验之常微分方程的数值解法（Python实现）

应用Python求解上篇函数在一系列离散点处的近似值。

# -*- coding: utf-8 -*-

import numpy as np
import random as r
import matplotlib.pyplot as plt
import math

p0=2
lam=5
#dp/dt=lam*p,p(0)=p0,求p(t)
x=[]
y=[]
for i in np.arange(0, 1.5, 0.001):
    x.append(i)
    y.append(p0*math.exp(lam*i))

plt.plot(x,y,linewidth=5)


#欧拉求解,向前分成300点分
xx=[0]
yy=[p0]
ds=1*1.0/100
for i in np.arange(0+ds, 1.5, ds):
    xx.append(i)
    yy.append(yy[-1]+lam*yy[-1]*ds)
    
plt.plot(xx,yy,'b',marker='+') 


#欧拉求解,向后分成300点分
xx2=[0]
yy2=[p0]
ds=1*1.0/100
for i in np.arange(0+ds, 1.5, ds):
    xx2.append(i)
    yy2.append(yy2[-1]*1.0/(1-lam*ds))
    
plt.plot(xx2,yy2,'r',marker='*') 

#欧拉求解,修正分成300点分
xx3=[0]
yy3=[p0]
ds=1*1.0/100
for i in np.arange(0+ds, 1.5, ds):
    xx3.append(i)
    yy3.append((yy3[-1]+ds/2*lam*yy3[-1])/(1-ds/2*lam))
    
plt.plot(xx3,yy3,'--','r',linewidth=2) 

plt.show()

运行结果：