实验4常微分方程数值解.ppt

《实验4常微分方程数值解.ppt》由会员分享，可在线阅读，更多相关《实验4常微分方程数值解.ppt（43页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、实验实验4 常微分方 程数值解凯 里 学 院 理 学 院主讲主讲:潘东云潘东云Experiments in Mathematics1为什么要学习微分方程数值解为什么要学习微分方程数值解 微分方程是研究函数变化规律的重要工具，有着广泛微分方程是研究函数变化规律的重要工具，有着广泛的应用。如的应用。如:物体的运动物体的运动,电路的电压电路的电压,人口增长的预测人口增长的预测 许多微分方程没有解析解，数值解法是求解的重要手许多微分方程没有解析解，数值解法是求解的重要手段，如段，如2实验实验4 4的基本内容的基本内容3.实际问题用微分方程建模，并求解实际问题用微分方程建模，并求解2.龙格龙格-库塔方法

2、的库塔方法的MATLAB实现实现*4.数值算法的收敛性、稳定性与刚性方程数值算法的收敛性、稳定性与刚性方程1.两个最常用的数值算法两个最常用的数值算法:欧拉（欧拉（Euler）方法）方法 龙格龙格-库塔（库塔（Runge-Kutta）方法）方法3实例实例1 1 海上缉私海上缉私 海防某部缉私艇上的雷达发现正东方向海防某部缉私艇上的雷达发现正东方向c海里处有一艘走私船正海里处有一艘走私船正以速度以速度a向正北方向行驶，缉私艇立即以最大速度向正北方向行驶，缉私艇立即以最大速度b(a)前往拦前往拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始终指截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始

3、终指向走私船。向走私船。建立任意时刻缉私艇位置及建立任意时刻缉私艇位置及 航线的数学模型航线的数学模型,并求解并求解;求出缉私艇追上走私船的时间。求出缉私艇追上走私船的时间。a北北bc艇艇船船4实例实例1 1 海上缉私海上缉私 建立坐标系如图建立坐标系如图:t=0 艇在艇在(0,0),船在船在(c,0);船速船速a,艇速艇速b 时刻时刻 t 艇位于艇位于P(x,y),船到达船到达 Q(c,at)模型模型:0yxcR(c,y)Q(c,at)P(x,y)b由方程无法得到由方程无法得到x(t),y(t)的解析解的解析解需要用数值解法求解需要用数值解法求解 5“常微分方程初值问题数值解常微分方程初值问

4、题数值解”的提法的提法不求解析解不求解析解 y=y(x)(无解析解或求解困难无解析解或求解困难)而在一系列离散点而在一系列离散点通常取等步长通常取等步长hy1y2yn6yP0 x0 x1x2x3 xy=y(x)y0P1P2P3欧 拉 方 法基本基本思路思路x取不同点取不同点各种欧拉公式各种欧拉公式向前欧拉公式向前欧拉公式显式公式显式公式7P3欧 拉 方 法向后欧拉公式向后欧拉公式 隐式公式隐式公式yP0 x0 x1x2x3 xy0y=y(x)P1P28向前欧拉公式向前欧拉公式向后欧拉公式向后欧拉公式二者平均得到二者平均得到梯形公式梯形公式仍为隐式公式，需迭代求解仍为隐式公式，需迭代求解将梯形公

5、式的迭代过程简化为两步将梯形公式的迭代过程简化为两步预测预测校正校正改进欧拉公式改进欧拉公式9龙格龙格-库塔方法库塔方法 向前，向后欧拉公式：向前，向后欧拉公式：用用xn,xn+1内个点的导内个点的导 数代替数代替 f(x,y(x)梯形公式，改进欧拉公式：梯形公式，改进欧拉公式：用用xn,xn+1内个点导数的平均值代替内个点导数的平均值代替 f(x,y(x)龙格龙格-库塔方法的基本思想库塔方法的基本思想在在xn,xn+1内多取几个点，将它们的导数加权内多取几个点，将它们的导数加权平均代替平均代替 f(x,y(x)，设法构造出设法构造出精度更高精度更高的计算公式。的计算公式。10常用的常用的(经

6、典经典)龙格龙格库塔库塔公式公式不足不足:收敛速度较慢收敛速度较慢L级?阶龙格龙格-库塔库塔方法的方法的一般形式一般形式使精度尽量高使精度尽量高4级4阶11微分方程组和高阶方程初值问题的数值解微分方程组和高阶方程初值问题的数值解 欧拉方法和龙格欧拉方法和龙格-库塔方法可直接推广到微分方程组库塔方法可直接推广到微分方程组向前欧拉公式向前欧拉公式 高阶方程需要先降阶为一阶微分方程组高阶方程需要先降阶为一阶微分方程组 12龙格龙格库塔方法的库塔方法的 MATLAB 实现实现 t,x=ode23(f,ts,x0,opt)3级级2阶龙格阶龙格-库塔库塔公式公式 t,x=ode45(f,ts,x0,opt

7、)5级级4阶龙格阶龙格-库塔库塔公式公式 f是待解方程写成是待解方程写成的函数的函数m文件：文件：function dx=f(t,x)dx=f1;f2;fn;ts=t0,t1,tf输出指定时刻输出指定时刻 t0,t1,tf 的函数值的函数值ts=t0:k:tf输出输出 t0,tf 内等分点处的函数内等分点处的函数值值x0为函数初值为函数初值(n维维)输出输出t=ts,x为相应函数值为相应函数值(n维维)optopt为选项，缺省时精度为：相对误差为选项，缺省时精度为：相对误差1010-3-3，绝对误差绝对误差1010-6-6,计算步长按精度要求自动调整计算步长按精度要求自动调整13微分方程的符号

8、解微分方程的符号解 求微分方程（组）的解析解(符号解）命令:dsolve(方程方程1,方程方程2,方程方程n,初始条件初始条件,自变量自变量)注意:在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分.任何D后所跟的字母为因变量，自变量可以指定或由系统规则选定为确省.例如，微分方程 应表达为：D2y=0.结 果：u=tan(t-c)14 解解 输入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)结 果 为:y=3e-2xsin（5x）15解解 输入命令：x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,

9、Dz=4*x-4*y+2*z,t)；x=simple(x)%将x化简 y=simple(y)z=simple(z)结 果 为：x=(c1-c2+c3+c2e-3t-c3e-3t)e2t y=-c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z=(-c1e-4t+c2e-4t+c1-c2+c3)e2t 16解：用解：用MATLABMATLAB编程如下：编程如下：function dx=fun113(t,x)dx=t+x;ts=0:0.1:1;x0=1;t,x=ode45(fun113,ts,x0);y=2*exp(t)-t-1;t,x,y,x-yplot(t,x

10、,r-,t,y,b-.),grid,gtext(x(t),gtext(y(t),title(微分方程数值解例微分方程数值解例);xlabel(time t);ylabel(solution x);legend(x,y);微分方程的数值解微分方程的数值解17微分方程的数值解微分方程的数值解解解:令 y1=x，y2=y11、建立m-文件fun114.m如下：function dx=fun114(t,x)dx=zeros(2,1);dx(1)=x(2);dx(2)=-x(2)/t+(1-0.25/t2)*x(1);2、取初值，输入命令：t,x=ode45(fun114,pi/2 2,pi/2-2/p

11、i plot(t,x(:,1),-)18实例实例1海上缉私海上缉私(续续)模型的数值解模型的数值解0yxc(x(t),y(t)ab设设:船速船速a=20(海里海里/小时小时)艇速艇速b=40(海里海里/小时小时)距离距离c=15(海里海里)求求:缉私艇的位置缉私艇的位置x(t),y(t)缉私艇的航线缉私艇的航线y(x)jisi.m,seajisi.m19 实例实例1海上缉私海上缉私(续续)模型的数值解模型的数值解 a=20,b=40,c=15走私船的位置走私船的位置x1(t)=c=15y1(t)=at=20t t=0.5时缉私艇追上走私船时缉私艇追上走私船 缉私艇的航线缉私艇的航线y(x)tx

12、(t)y(t)0000.051.99840.06980.103.98540.29240.155.94450.69060.207.85151.28990.259.67052.11780.3011.34963.20050.3512.81704.55520.4013.98066.17730.4514.74518.02730.5015.00469.9979y1(t)01.02.03.04.05.06.07.08.09.010.020实例实例1 海上缉私海上缉私(续续)模型的数值解模型的数值解 设设b，c不变，不变，a变大变大为为30,35,接近接近40,观察解的变化观察解的变化:a=35,b=40,c

13、=15t=?缉私艇追上走私船缉私艇追上走私船 t x(t)y(t)y1(t)0 0 0 00.13.95610.50583.50.27.59282.13087.00.310.52404.828310.50.412.53848.275514.00.513.755112.083017.51.214.998640.016442.01.314.999644.016545.51.415.011748.018349.01.515.002352.014652.51.614.986655.948656.0累积误差较大累积误差较大提高精度提高精度!21实例实例1海上缉私海上缉私(续续)模型的数值解模型的数值解

14、a=35,b=40,c=15opt=odeset(RelTol,1e-6,AbsTol,1e-9);t,x=ode45(jisi,ts,x0,opt);t=1.6时缉私艇追上走私船时缉私艇追上走私船 缉私艇的航线缉私艇的航线y(x)判断判断“追上追上”的有效方的有效方法法?t x(t)y(t)y1(t)0 0 0 00.13.9561040.5058133.50.27.5928222.1306787.00.310.5219214.82930810.50.412.5394548.26984014.00.513.75397412.07534417.51.214.99961640.00000542.

15、01.314.99996344.00000545.51.414.99999348.00000549.01.514.99999852.00000552.51.615.00002055.99993156.022实例实例1 海上缉私海上缉私(续续)模型的解析解模型的解析解 23实例实例1海上缉私海上缉私(续续)模型的解析解模型的解析解 缉私艇的航线缉私艇的航线y y(x x)的解析解的解析解x=c时时 缉私艇追上走私船的缉私艇追上走私船的y y坐标坐标 缉私艇追上走私船的时间缉私艇追上走私船的时间:a=20,b=40,c=15 t1=0.5 a=35,b=40,c=15 t1=1.6 24微分方程数

16、值解法的误差分析微分方程数值解法的误差分析数值解法数值解法:计算微分方程精确解计算微分方程精确解y(xn)的近似值的近似值yn按照步长按照步长h一步步计算一步步计算,每步都有误差每步都有误差;每一步的误差会逐步积累每一步的误差会逐步积累,称称累积误差累积误差.讨论计算一步出现的误差讨论计算一步出现的误差假定假定yn=y(xn),估计估计yn+1的误差的误差:y(xn+1)-yn+1局部截断误差局部截断误差 25误差分析误差分析估计欧拉公式的局部截断误差估计欧拉公式的局部截断误差 y(xn+1)在在xn处作处作Taylor展开展开:向前欧拉公式向前欧拉公式yn=y(xn)局部截断误差主项为局部截

17、断误差主项为26误差分析误差分析估计欧拉公式的局部截断误差估计欧拉公式的局部截断误差 向后欧拉公式向后欧拉公式局部截断误差主项为局部截断误差主项为梯形梯形公式公式向前、向后欧拉公式的平均向前、向后欧拉公式的平均 xn xn+1 x y Pn A Q B27误误差差分分析析算法算法精度的阶精度的阶的定义的定义 一个算法的局部截断误差为一个算法的局部截断误差为O(hp+1)该算法具有该算法具有p阶精度阶精度 局部截断误差局部截断误差精度精度向前欧拉公式向前欧拉公式O(h2)1阶阶向后欧拉公式向后欧拉公式O(h2)1阶阶梯形公式梯形公式O(h3)2阶阶改进欧拉公式改进欧拉公式O(h3)2阶阶经典龙格

18、经典龙格-库塔公式库塔公式O(h5)4阶阶28实实 例例 2弱弱 肉肉 强强 食食 问问题题自然界中同一环境下两个种群之间的生存方式自然界中同一环境下两个种群之间的生存方式 相互竞争相互竞争 相互依存相互依存 弱肉强食弱肉强食 弱肉弱肉强食强食种群甲靠丰富的自然资源生存种群甲靠丰富的自然资源生存 食饵食饵(Prey)(Prey)种群乙靠捕食种群甲为生种群乙靠捕食种群甲为生 捕食者捕食者(Predator)(Predator)两个种群的数量如何演变两个种群的数量如何演变?历史背景历史背景一次世界大战期间地中海渔业一次世界大战期间地中海渔业的捕捞量下降的捕捞量下降(食用鱼和鲨鱼同时捕捞食用鱼和鲨鱼

19、同时捕捞)，但，但是其中是其中鲨鱼的比例却增加，为什么？鲨鱼的比例却增加，为什么？29实实 例例 2弱弱 肉肉 强强 食食 模型模型食饵食饵(甲甲)的密度的密度x(t),捕食者捕食者(乙乙)的密度的密度y(t)甲独立生存的增长率甲独立生存的增长率r r乙使甲的增长率减小乙使甲的增长率减小,减小量与减小量与 y y 成正比成正比乙独立生存的死亡率乙独立生存的死亡率d甲使乙的死亡率减小甲使乙的死亡率减小,减小量与减小量与 x成正比成正比Volterra模型模型 x(t),y(t)无解析解无解析解30实例实例2弱肉强食弱肉强食 模型的数值解模型的数值解 猜测猜测 x(t),y(t)是周期函数是周期函

20、数;y(x)是封闭曲线是封闭曲线 数值积分计算一个周期的平均值：数值积分计算一个周期的平均值：shier.m,shier1.m31实例实例2弱肉强食弱肉强食 模型的解析解模型的解析解 c由初始条件确定由初始条件确定 相轨线是封闭曲线相轨线是封闭曲线(c在一定范围内在一定范围内)求求x(t),y(t)一周期的平均值一周期的平均值:可以可以证明证明x(t),y(t)是周期函数是周期函数(周期记作周期记作T)32实例实例2弱肉强食弱肉强食 模型的解析解模型的解析解 x(t),y(t)一周期的平均值一周期的平均值:r 食饵增长率食饵增长率a 捕食者对食饵的捕获能力捕食者对食饵的捕获能力 d 捕食者死亡

21、率捕食者死亡率b 食饵对捕食者的喂养能力食饵对捕食者的喂养能力 结果解释结果解释与计算结果同与计算结果同既相互制约既相互制约又相互依存又相互依存33T2T3T4T1PT1 T2 T3 T4x(t)的的“相位相位”领先领先 y(t)进一步分析进一步分析初值初值相轨线的方向相轨线的方向34一次大战期间地中海渔业的捕捞量下降，一次大战期间地中海渔业的捕捞量下降，但是其中但是其中鲨鱼的比例却在增加，为什么？鲨鱼的比例却在增加，为什么？rr-1,dd+1捕捞捕捞战时战时捕捞捕捞rr-2,dd+2,2 0微分方程稳定微分方程稳定 向前欧拉公式向前欧拉公式向后欧拉公式向后欧拉公式经典经典龙格龙格-库塔公式库

22、塔公式算法稳定算法稳定(特征根特征根-)h任意任意39 刚性现象与刚性方程刚性现象与刚性方程振动系统或电路系统的数学模型振动系统或电路系统的数学模型 现象现象 k=2000.5,r=1000,a=1,b=-1999.5,f(t)=1 瞬态解与稳态解瞬态解与稳态解 e-2000t 快瞬态解快瞬态解 e-t/2 慢瞬态解慢瞬态解 计算到计算到t=0.005时已衰减到时已衰减到4.5 10-5 计算到计算到t=20时才衰减到时才衰减到4.5 10-5 精度达到精度达到10-4需算到需算到t=20(由慢瞬态解由慢瞬态解=1/2决定决定)选取步长选取步长h由快瞬态解由快瞬态解=2000决定决定h2.78

23、5/2000=0.0014 龙格龙格-库塔公式库塔公式t=20需需14286步步 快快、慢瞬态解的特征根相差悬殊慢瞬态解的特征根相差悬殊 刚性现象刚性现象(Stiff)求稳求稳态解态解40 刚性现象与刚性方程刚性现象与刚性方程刚性刚性方程方程 振动、电路及化学反应中的线性常系数方程组振动、电路及化学反应中的线性常系数方程组 A的特征根的特征根 k(k=1,2,n)的实部的实部Re(k)10刚性非线性方程组线性常系数方程组传统的数值方法无能为力传统的数值方法无能为力 线性化方法刚性方程刚性方程 41 刚性现象与刚性方程刚性现象与刚性方程刚性方程的刚性方程的MATLAB求解求解 ode23,ode45 解刚性方程的困难解刚性方程的困难步长自动变小步长自动变小 计算时间很长计算时间很长 求解刚性方程的命令求解刚性方程的命令:ode23s,ode15s 等等(用法相同用法相同)例例特征根特征根 1=-1,2=-106刚性比刚性比 s=106 解析解解析解Stiff.m,stiff1.m42布置实验布置实验目的目的1.用用MATLAB软件掌握求微分方程数值解的方法，软件掌握求微分方程数值解的方法，并对结果作初步分析；并对结果作初步分析；2.通过实例学习用微分方程模型解决简化的实际问通过实例学习用微分方程模型解决简化的实际问题。题。内容内容1，5，7，843