重庆工商大学数学模型与数学实验课件第08讲 常微分方程数值解

大学数学实验ExperimentsinMathematics实验4常微分方程数值解《数学建模与数学实验》---李焕荣为什么要学习微分方程数值解•微分方程是研究函数变化规律的重要工具，有着广泛的应用。如:物体的运动,电路的电压,人口增长的预测•许多微分方程没有解析解，数值解法是求解的重要手段，如2,dyyxdx()()xtaxyytaxyby实验4的基本内容3.实际问题用微分方程建模，并求解2.龙格-库塔方法的MATLAB实现*4.数值算法的收敛性、稳定性与刚性方程1.两个最常用的数值算法:•欧拉（Euler）方法•龙格-库塔（Runge-Kutta）方法实例1海上缉私海防某部缉私艇上的雷达发现正东方向c海里处有一艘走私船正以速度a向正北方向行驶，缉私艇立即以最大速度b(a)前往拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始终指向走私船。•建立任意时刻缉私艇位置及航线的数学模型,并求解;•求出缉私艇追上走私船的时间。a北bc艇船实例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)bcos,sindxdybbdtdt2222()()()()()()dxbcxdtcxatydybatydtcxaty由方程无法得到x(t),y(t)的解析解需要用数值解法求解“常微分方程初值问题数值解”的提法00'(,),()yfxyyxy设的解y=y(x)存在且唯一不求解析解y=y(x)(无解析解或求解困难)12nxxx而在一系列离散点()(1,2,)nyxnn求的近似值y通常取等步长h0nxxnh0x0y)(xyyyxy1y2yn)(nxy1x2xnxyP0x0x1x2x3xy=y(x)y0P1P2P3欧拉方法00'(,),()yfxyyxy基本思路1[,]nnxx在小区间上1(,)[,]nnfxyxx中的x取内的某一点1'[()()]/,nnyyxyxh11()()(,()),[,]nnnnyxyxhfxyxxxxx取不同点各种欧拉公式nxx取左端点1()()(,())nnnnyxyxhfxyx1(,),0,1,nnnnyyhfxyn向前欧拉公式显式公式11(),()nnnnyyxyyxP3欧拉方法11()()(,()),[,]nnnnyxyxhfxyxxxx11(),()nnnnxyyxyyx取右端点，111(,),0,1,nnnnyyhfxyn向后欧拉公式隐式公式yP0x0x1x2x3xy0y=y(x)P1P2n+1y右端未知，需迭代求解(0)1(,)nnnnyyhfxy初值(1)()111(,)0,1,2,,0,1,2,kknnnnyyhfxykn()11limknnkyy向前欧拉公式向后欧拉公式1(,)nnnnyyhfxy111(,)nnnnyyhfxy二者平均得到梯形公式111[(,)(,)],0,1,2nnnnnnhyyfxyfxyn仍为隐式公式，需迭代求解将梯形公式的迭代过程简化为两步1(,)nnnnyyhfxy预测111[(,)(,)],0,1,2nnnnnnhyyfxyfxyn校正改进欧拉公式龙格-库塔方法11()()(,()),[,]nnnnyxyxhfxyxxxx•向前，向后欧拉公式：用[xn,xn+1]内１个点的导数代替f(x,y(x))•梯形公式，改进欧拉公式：用[xn,xn+1]内２个点导数的平均值代替f(x,y(x))龙格-库塔方法的基本思想在[xn,xn+1]内多取几个点，将它们的导数加权平均代替f(x,y(x))，设法构造出精度更高的计算公式。常用的(经典)龙格—库塔公式不足:收敛速度较慢111222111(,)(,)(,),3,4,,LnniiinnnniininiijjjyyhkkfxykfxchychkkfxchychakiLL级?阶龙格-库塔方法的一般形式112341213243(22)6(,)(/2,/2)(/2,/2)(,)nnnnnnnnnnhyykkkkkfxykfxhyhkkfxhyhkkfxhyhk1,10,1,,111ijijiLiiijiiacac满足使精度尽量高4级4阶微分方程组和高阶方程初值问题的数值解欧拉方法和龙格-库塔方法可直接推广到微分方程组0000(,,)(,,)(),()yfxyzzgxyzyxyzxz向前欧拉公式),,(yyxfy1yy),,(21221yyxfyyy高阶方程需要先降阶为一阶微分方程组,2,1,0),,(),,(11nzyxhgzzzyxhfyynnnnnnnnnn龙格—库塔方法的MATLAB实现TnTnfffxxxxtxxtftx),(,),(,)(),,()(1100[t,x]=ode23(@f,ts,x0)3级2阶龙格-库塔公式[t,x]=ode45(@f,ts,x0)5级4阶龙格-库塔公式f是待解方程写成的函数m文件：functiondx=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维)缺省精度(相对误差10-3，绝对误差10-6),计算步长按精度要求自动调整.实例1海上缉私(续)模型的数值解2222()()()()()()dxbcxdtcxatydybatydtcxaty(0)0,(0)0xy0yxc(x(t),y(t))ab设:船速a=20(海里/小时)艇速b=40(海里/小时)距离c=15(海里)求:缉私艇的位置x(t),y(t)缉私艇的航线y(x)实例1海上缉私(续)模型的数值解a=20,b=40,c=15走私船的位置x1(t)=c=15y1(t)=at=20tt=0.5时缉私艇追上走私船051015200246810xy缉私艇的航线y(x)tx(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.0刚性现象与刚性方程bxaxrktfrxxkx)0(,)0(),0,()(振动系统或电路系统的数学模型现象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.785/2000=0.0014龙格-库塔公式t=20需14286步快、慢瞬态解的特征根相差悬殊2000/2()1ttxtee刚性现象(Stiff)求稳态解刚性现象与刚性方程刚性方程的MATLAB求解ode23,ode45解刚性方程的困难步长自动变小计算时间很长求解刚性方程的命令:ode23s,ode15s等(用法相同)1126621266122(101)(102)(0)10/4,(0)10/41/2xxxxxxxx例特征根1=-1,2=-106刚性比s=1066661016610210()(1)410(101)()(1)42ttttxteextee解析解Stiff.m,stiff1.m布置实验目的1.用MATLAB软件掌握求微分方程数值解的方法，并对结果作初步分析；2.通过实例学习用微分方程模型解决简化的实际问题。内容5