数值分析方法 课件 第六章 常微分方程数值解法

数值分析方法面向“四新”人才培养普通高等教育系列教材第六章常微分方程数值解法目录/Contents6.1

认识微分方程

6.2

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

6.3

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

6.4

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

6.5

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

目录/Contents6.6

线性多步方法

6.7

微分方程组的刚性问题

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

6.9

微分方程计算机实验

含未知函数及其导数的方程叫做常微分方程如果未知函数是多元函数且方程中含有未知函数的偏导数则称其为偏微分方程常微分方程中未知函数导数的阶数称为微分方程的阶数含有任意常数的函数满足微分方程且任意常数的个数与微分方程的阶数相等称为该常微分方程的通解既满足微分方程且适合初值条件的解称为该微分方程初值问题的特解6.1

认识微分方程6.1.1微分方程模型举例

因此求解满足

解的问题称为微分方程的初值问题。例6.1.2振动模型

介质中质量为m的质点，假定处在弹性约束之下作一维振动(即仅需一个位置参数就可完全描述质点状态的运动)，通常以弹簧作为这类一维弹性振动的代表模型（图6-1）。已知质点在介质中运动所受阻力与质点速度成正比，根据牛顿第二定律，

SIS模型只把人群分为S类和I类，患者病愈后并未产生免疫力，仍未S类成员。SIS模型中各类人员的相互作用关系如图6-2所示。NSNI

SIR模型考虑S类、I类和R类人群，且I类人病愈后具有免疫能力，不会被传染而进入I类。S，I，R三类人口的关系如图6-3所示。NSNINR

方程组（6.1.5）常称为SIR模型。例6.1.4洛伦兹方程

考虑底部受热的均匀厚度的水平流体层。设底部温度高于顶部温度，且流体受到重力作用。当底部流体受热膨胀时，下面流体密度小于上面流体密度，于是在重力场的作用下，上面流体有下沉的趋势，而下面流体有上浮的趋势。当温度梯度较小时，由于流体粘性的阻碍，流体仍处于静止状态，在流体内部只有热传导，温度随高度呈线性变化。然而当温度梯度较大时，流体的静止平衡状态失稳，出现对流现象。这一流体的运动可以用一个复杂的偏微分方程表示。

6.1.2微分方程数值解一般地，对于一阶常微分方程的初值问题

（6.1.8）对于一个由多个变量构成微分方程系统，或称微分方程组，本质上解法是一样的。一般地，一个由n个方程组成的方程组，具有如下形式

谢谢数值分析方法主编

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

认识微分方程

6.2

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

6.3

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

6.4

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

6.5

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

目录/Contents6.6

线性多步方法

6.7

微分方程组的刚性问题

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

6.9

微分方程计算机实验

6.2微分方程初值问题的Euler方法6.2.1显式Euler公式

xEuler法y精确解01.000001.000000.11.100001.089960.21.180941.15940.31.242271.207070.41.282771.230960.51.300611.228220.61.293251.194970.71.257361.125890.81.188541.013400.91.080820.8457811.00.9258020.6016496.2.2隐式Euler公式与改进Euler公式

例6.2.2

用改进Euler公式（6.2.10）求解初值问题（6.2.3）。xEuler法y精确解改进的Euler法y01.000001.0000010.11.100001.089961.090470.21.180941.15941.160470.31.242271.207071.20880.41.282771.230961.233480.51.300611.228221.231740.61.293251.194971.199760.71.257361.125891.132360.81.188541.013401.022160.91.080820.8457810.8578951.00.9258020.6016490.619366谢谢数值分析方法主编

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

认识微分方程

6.2

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

6.3

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

6.4

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

6.5

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

目录/Contents6.6

线性多步方法

6.7

微分方程组的刚性问题

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

6.9

微分方程计算机实验

6.3初值问题数值解的误差与稳定性分析6.3.1误差分析

（6.3.1）

6.3.2收敛性与稳定性分析1.数值方法的收敛性

2.数值方法的稳定性

谢谢数值分析方法主编

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

认识微分方程

6.2

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

6.3

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

6.4

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

6.5

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

目录/Contents6.6

线性多步方法

6.7

微分方程组的刚性问题

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

6.9

微分方程计算机实验

6.4.微分方程初值问题的Runge-Kutta法6.4.1Runge-Kutta法的基本思想与二阶Runge-Kutta法

例6.4.1

用二阶Runge-Kutta公式（6.4.2）求解初值问题（6.2.3）。

二阶Runge-Kutta公式（6.4.2）所求解与改进Euler方法所求解具有基本一致的误差，图6-4（右）为两种算法与精确解比较的结果，其中菱形与圆点数据分别表示相同节点处精确解分别与二阶Runge-Kutta公式（6.4.2）所求解和改进Euler方法所求解之差。可以看出二种算法有微小差别。6.4.2三、四阶Runge-Kutta法

计算的结果见图6-5，其中实心方块为精确解，实心三角为三阶Runge-Kutta公式，空心三角为二阶Runge-Kutta公式所求结果，空心方块为改进Euler方法所求解。

6.4.3隐式Runge-Kutta法

谢谢数值分析方法主编

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

认识微分方程

6.2

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

6.3

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

6.4

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

6.5

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

目录/Contents6.6

线性多步方法

6.7

微分方程组的刚性问题

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

6.9

微分方程计算机实验

6.5非线性微分方程组初值问题的

龙格-库塔法

谢谢数值分析方法主编

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

认识微分方程

6.2

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

6.3

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

6.4

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

6.5

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

目录/Contents6.6

线性多步方法

6.7

微分方程组的刚性问题

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

6.9

微分方程计算机实验

6.6线性多步方法

6.6.1线性多步方法的构造

（6.6.1）

6.6.2线性多步方法的应用及预测-校正方法1.线性多步方法的应用

2.预测-校正方法

2.线性多步方法的收敛性与稳定性

谢谢数值分析方法主编

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

认识微分方程

6.2

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

6.3

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

6.4

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

6.5

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

目录/Contents6.6

线性多步方法

6.7

微分方程组的刚性问题

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

6.9

微分方程计算机实验

6.7微分方程组的刚性问题

求解微分方程组时，经常出现解的分量数量级差别很大的情形，这给数值求解带来了较大困难，这种问题称为刚性问题。

步数nxyz10202-7.0354920.1041-6.828563-261.236134014.-8989.144-7.13065*10^147.47272*10^26-3.39266*10^255-4.98155*10^1044.42368*10^172-2.00838*10^171

谢谢数值分析方法主编

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

认识微分方程

6.2

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

6.3

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

6.4

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

6.5

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

目录/Contents6.6

线性多步方法

6.7

微分方程组的刚性问题

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

6.9

微分方程计算机实验

6.8二阶微分方程的边值问题二阶微分方程也可以看作是微分方程组，如（6.1.2）可以化作（6.1.10），可以利用微分方程组的数值方法求解实际应用中较为常见的微分方程的初值问题，但对二阶微分方程而言，所谓的边值问题在科学技术问题中也常常出现，其求解问题与上述初值问题不同，本节简单讨论两种本方法。6.8.1二阶微分方程边值问题的打靶法

6.8.2二阶线性微分方程边值问题的差分法

固体表面扩散层液体生物膜

谢谢数值分析方法主编

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