数值分析方法 课件 6.9 微分方程计算机实验

数值分析方法主编

李冬果李林高磊首都医科大学生物医学工程学院智能医学工程学学系面向“四新”人才培养普通高等教育系列教材第六章常微分方程数值解法目录/Contents6.1

认识微分方程

6.2

微分方程初值问题的Euler方法

6.3

微分方程初值问题数值解的误差与稳定性分析

6.4

微分方程初值问题的Runge-Kutta法

6.5

非线性微分方程组初值问题的龙格-库塔法

目录/Contents6.6

线性多步方法

6.7

微分方程组的刚性问题

6.8二阶微分方程的边值问题

6.9

微分方程计算机实验

6.9微分方程计算机实验6.9.1显式Euler公式和改进Euler公式的实现我们可以通过以下代码实现Euler法：defeuler(f,x0,y0,b,h):

y=[y0]

x=[x0]

while(x0+h)<=b:

y0=h*f(x0,y0)+y0

x0=x0+h

y.append(y0)

x.append(x0)

returny,x其中输入参数:函数f即为所求的初值问题中的导函数，x0为初值点，y0为初值点处对应的函数值，b为积分区间的终点，而h则是步长。返回值为两个列表，其中分别存储计算出的微分方程解的函数值，和与该函数值对应的自变量值。类似于Euler公式，将递推公式进行修改，就得到了改进Euler法的求解函数defeuler_improved(f,x0,y0,b,h):

y=[y0]

x=[x0]

while(x0+h)<=b:

x1=x0+h

yp=y0+h*f(x0,y0)

yc=y0+h*f(x1,yp)

y0=(yp+yc)/2

x0=x1

y.append(y0)

x.append(x0)

returny,x其参数与之前Euler法的定义相同。例6.9.1用Euler公式（6.2.2）和改进Euler公式（6.2.10）求解初值问题（6.2.3）。解：

为了实现指数运算，需要导入其他工具包，这里使用numpy工具包。importnumpyasnp使用lambda定义导函数fun_6_2_1和解函数real_ans_6_2_1fun_6_2_1=lambdax,y:y-3*x/y**(1/3)

real_ans_6_2_1=lambdax:(9+12*x-5*np.exp(4*x/3))**(3/4)/2**(3/2)应用之前定义的函数实现Euler公式和改进Euler公式，将函数值分别返回给y1和y2[y1,x1]=euler(fun_6_2_1,0,1,1,0.1)

[y2,x2]=euler_improved(fun_6_2_1,0,1,1,0.1)计算真实的函数值y_real：

y_real=[real_ans_6_2_1(x)forxinx1]

打印结果：print("x\tEuler法y\t真实值\t改进Euler法y")

forkinrange(11):

print(f"{x1[k]:.1f}\t{y1[k]:.6f}\t{y_real[k]:.6f}\t{y2[k]:.6f}")结果见表6-7。6.9.2四阶Runge-Kutta法的实现类似于之前Euler和改进Euler法的实现，可以根据公式（6.4.7）实现四阶的Runge-Kutta法，可以如下定义函数defrk4(f,x0,y0,b,h):

x=[x0]

y=[y0]

while(x0+h)<=b:

x1=x0+h

k1=f(x0,y0)

k2=f(x0+h/2,y0+h/2*k1)

k3=f(x0+h/2,y0+h/2*k2)

k4=f(x1,y0+h*k3)

y0=y0+h*(k1+2*k2+2*k3+k4)/6

x0=x1

y.append(y0)

x.append(x0)

returny,x其参数与之前Euler法的定义相同。

可以使用matplotlib工具包来绘制函数图像，为了方便数据处理，同时引入numpy工具包importnumpyasnp

importmatplotlib.pyplotasplt

x1_array=np.array(x1)

Y_array=np.array(Y)

markers=["s","o","s","^"]

markerfillstyles=["none","full","full","full"]

plt.figure(figsize=(10,6))

forkinrange(4):

plt.plot(x1_array[:200:5],Y_array[k,:200:5],

marker=markers[k],fillstyle=markerfillstyles[k],

linestyle='',color='blue')

plt.axhline(5,0,10,linestyle='-')

plt.axhline(2,0,10,linestyle='--')

plt.show()结果见图6-6。6.9.3方程组的四阶Runge-Kutta法实现根据公式（6.5.2），可以得到方程组的数值计算方法，可以编写Python函数如下。求解由2个方程组成的方程组的函数rk4_2:defrk4_2(f,g,t0,x0,y0,b,h):

t=[t0]

x=[x0]

y=[y0]

while(t0+h)<=b:

t1=t0+h

k1=f(t0,x0,y0)

l1=g(t0,x0,y0)

k2=f(t0+h/2,x0+h/2*k1,y0+h/2*l1)

l2=g(t0+h/2,x0+h/2*k1,y0+h/2*l1)

k3=f(t0+h/2,x0+h/2*k2,y0+h/2*l2)

l3=g(t0+h/2,x0+h/2*k2,y0+h/2*l2)

k4=f(t1,x0+h*k3,y0+h*l3)

l4=g(t1,x0+h*k3,y0+h*l3)

x0=x0+h*(k1+2*k2+2*k3+k4)/6

y0=y0+h*(l1+2*l2+2*l3+l4)/6

t0=t1

x.append(x0)

y.append(y0)

t.append(t0)

returnx,y,t

同理可以计算无阻尼假设下的运动轨迹：g2=lambdat,x,y:-1.2*x-0.1*y

x2,y2,t2=rk4_2(f,g2,t0,x0,y0,b,h)

可以使用matplotlib绘制图像：importmatplotlib.pyplotasplt

plt.figure(figsize=(14,6))

plt.subplot(1,2,2)

plt.plot(t1,x1)

plt.subplot(1,2,1)

plt.plot(t2,x2)

plt.show()由此得到无阻尼和正弦阻尼情形下的震动图像（见图6-7）。

相应Python程序如下：

si,b0,r=10,8/3,28

f=lambdat,x,y,z:si*(y-x)

g=lambdat,x,y,z:r*x-y-x*z

p=lambdat,x,y,z:y*x-b0*z

t0,x0,y0,z0,b,h=0,0,2,0,34,0.01

x1,y1,z1,t1=rk4_3(f,g,p,t0,x0,y0,z0,b,h)修改参数：y0=-2

x2,y2,z2,t2=rk4_3(f,g,p,t0,x0,y0,z0,b,h)