《MATLAB及其在理工课程中的应用常微分方程4》由会员分享，可在线阅读，更多相关《MATLAB及其在理工课程中的应用常微分方程4（17页珍藏版）》请在金锄头文库上搜索。
1、 常微分方程常微分方程主要内容主要内容学习学习MATLABMATLAB有关求解常微分方程（组）的指令；有关求解常微分方程（组）的指令；刚性方程组问题；刚性方程组问题；实际问题的微分方程建模及求解方法；实际问题的微分方程建模及求解方法；解常微分方程主要解常微分方程主要MATLABMATLAB指令指令pode45 ode45 四、五阶四、五阶Runge-kuttaRunge-kutta法法pode23 ode23 二、三阶二、三阶Runge-kuttaRunge-kutta法法pode113 ode113 多步多步AdamsAdams算法算法podesetodeset 解解odeode选项设置选项
2、设置pode23t ode23t 适度刚性问题梯形算法适度刚性问题梯形算法pode15s ode15s 刚性方程组多步刚性方程组多步GearGear法法pode23s ode23s 刚性方程组二阶刚性方程组二阶RosebrockRosebrock法法pode23tb ode23tb 刚性方程组低精度算法刚性方程组低精度算法pbvpinitbvpinit 边值问题的预估解边值问题的预估解pbvp4c bvp4c 边值问题解法边值问题解法pdevaldeval 微分方程解的求值微分方程解的求值预备知识：常微分方程预备知识：常微分方程n微分方程的概念微分方程的概念常微分方程常微分方程偏微分方程偏微分
3、方程线性常微分方程线性常微分方程n初等积分法初等积分法( (分离变量法、积分因子法、常数变易法、降阶法等分离变量法、积分因子法、常数变易法、降阶法等) )n常系数线性微分方程常系数线性微分方程一阶变系数常微分方程：叠加原理求相应齐次微分方程解一阶变系数常微分方程：叠加原理求相应齐次微分方程解+ +一一个特解个特解高阶常系数微分方程：求相应齐次微分方程基本解高阶常系数微分方程：求相应齐次微分方程基本解+ +常数变易常数变易法求特解法求特解n初值问题数值解初值问题数值解注：高阶常微分方程初值问题需先转化为一阶常微分方程组注：高阶常微分方程初值问题需先转化为一阶常微分方程组一、初值问题求解一、初值问
4、题求解n常用调用格式：常用调用格式： t,yt,y=ode45(odefun,tspan,y0)=ode45(odefun,tspan,y0)u 参数说明：参数说明：uodefunodefun 表示表示f(t,yf(t,y) )的的函数句柄或函数句柄或inlineinline函数函数，t t是标量，是标量，y y是是标量或向量；标量或向量；utspantspan 如果是二维向量如果是二维向量t0,tft0,tf，表示自变量初值，表示自变量初值t0t0和和tftf；如；如果是高维向量果是高维向量t0,t1,t0,t1,tntn ，则表示输出结点列向量；，则表示输出结点列向量；uy0 y0 表示初
5、值向量表示初值向量y0y0；ut t 表示结点列向量表示结点列向量(t0,t1,(t0,t1,tn)Ttn)T; ;uy y 表示数值解矩阵，每一列对应表示数值解矩阵，每一列对应y y的一个分量；的一个分量；n完整调用格式：完整调用格式： t,yt,y=ode45(odefun,tspan,y0=ode45(odefun,tspan,y0，options,p1,p2,)options,p1,p2,)uoptions options 为计算参数设置（如精度要求等），默认为计算参数设置（如精度要求等），默认表示；表示；up1p1，p2p2， 附加传递参数，附加传递参数，odefunodefun表示
6、为表示为f(t,y,p1,p2,)f(t,y,p1,p2,)一、初值问题求解一、初值问题求解pode45ode45是最常用的求解微分方程的指令。采用变步长是最常用的求解微分方程的指令。采用变步长四、五阶龙格四、五阶龙格- -库塔法，适合高精度问题；库塔法，适合高精度问题；pode23ode23与与ode45ode45类似，但精度低一些；类似，但精度低一些；pode113ode113采用多步法，高低精度均适合；采用多步法，高低精度均适合；注：此三指令对于刚性方程组不宜采用。注：此三指令对于刚性方程组不宜采用。pode23tode23t，ode23s, ode23tb, ode15sode23s,
7、 ode23tb, ode15s都是求解刚性方都是求解刚性方程组的指令。程组的指令。一、初值问题求解一、初值问题求解例题：例题：1.1.解微分方程解微分方程odefunodefun=inline(y-2*=inline(y-2*t/y,t,yt/y,t,y););t,yt,y=ode45(odefun,0,4,1); =ode45(odefun,0,4,1); t,yt,yplot(t,y,oplot(t,y,o-)-)ans =ans =0 1.00000 1.00000.0502 1.04900.0502 1.04900.1005 1.09590.1005 1.09590.1507 1.1
8、4080.1507 1.1408.3.9502 2.98393.9502 2.98394.0000 3.00064.0000 3.0006ode45(odefun,0,4,1);ode45(odefun,0,4,1);t,yt,y=ode45(odefun,0:1:4,1);t,y=ode45(odefun,0:1:4,1);t,y一、初值问题求解一、初值问题求解例题：例题：2.2.解微分方程组解微分方程组解：将变量解：将变量x x，y y合写为向量变量合写为向量变量x x，做，做M M函数函数eg6_3fun.meg6_3fun.m%M%M函数函数eg6_3fun.meg6_3fun.mfu
9、nction f=eg6_3funfunction f=eg6_3fun（t t，x x）f(1)=-x(1)3-x(2);f(1)=-x(1)3-x(2);f(2)=x(1)-x(2)3;f(2)=x(1)-x(2)3;f=f(:);f=f(:);clear;t,xclear;t,x=ode45(=ode45(eg6_3fun,0 30,1;0.5);eg6_3fun,0 30,1;0.5);subplot(1,2,1);plot(t,x(:,1),t,x(:,2),:);subplot(1,2,1);plot(t,x(:,1),t,x(:,2),:);subplot(1,2,2);plot
10、(x(:,1), x(:,2);subplot(1,2,2);plot(x(:,1), x(:,2);一、初值问题求解一、初值问题求解例题：例题：3.3.解微分方程组（竖直加热板的自然对流）解微分方程组（竖直加热板的自然对流）已知：当已知：当=0=0时，时，解：首先引入辅助变量解：首先引入辅助变量将方程组化为一阶方程组将方程组化为一阶方程组一、初值问题求解一、初值问题求解先写先写M M函数函数eg6_4fun.meg6_4fun.m% M% M函数函数eg6_4fun.meg6_4fun.mfunction f=eg6_4fun(t,y)function f=eg6_4fun(t,y)f=y(
11、2);y(3);-3*y(1)*y(3)+2*y(2)2-y(4);f=y(2);y(3);-3*y(1)*y(3)+2*y(2)2-y(4); y(5);-2.1*y(1)*y(5); y(5);-2.1*y(1)*y(5);y0=0,0,0.68,1,-0.5;y0=0,0,0.68,1,-0.5;t,yt,y=ode45(eg6_4fun,0 5,y0);=ode45(eg6_4fun,0 5,y0);plot(t,y(:,1),t,y(:,4),:)plot(t,y(:,1),t,y(:,4),:)二、边值问题解法二、边值问题解法n常微分方程一阶方程组边值问题常微分方程一阶方程组边值问
12、题MATLABMATLAB标准形式：标准形式：n调用格式调用格式uSinitSinit= =bvpinit(tinit,yinitbvpinit(tinit,yinit) ) 由在粗略结点由在粗略结点tinittinit的预估解的预估解yinityinit生成粗略解网格生成粗略解网格sinitsinituSol=bvp4c(odefun,bcfun,sinit)Sol=bvp4c(odefun,bcfun,sinit) odefunodefun是微分方程组函数，是微分方程组函数，bcfunbcfun为边值条件函数，为边值条件函数，sol.xsol.x为求为求解结点，解结点，sol.ysol.y
13、是是y(ty(t) )的数值解的数值解uSxSx= =deval(sol,tideval(sol,ti) ) 计算由计算由bvp4cbvp4c得到的解在得到的解在titi的值的值二、边值问题解法二、边值问题解法例题：求解边值问题例题：求解边值问题解：首先改写成方程组解：首先改写成方程组 边界条件为边界条件为二、边值问题解法二、边值问题解法求解用求解用M M函数函数eg6_5fun.meg6_5fun.m% M% M函数函数eg6_5fun.meg6_5fun.mclear;closeclear;close; ;sinitsinit=bvpinit(0:4,1;0)=bvpinit(0:4,1;
14、0)odefunodefun=inline(y(2);-abs(y(1),t,y);=inline(y(2);-abs(y(1),t,y);bcfunbcfun=inline(ya(1);yb(1)+2,ya,yb);=inline(ya(1);yb(1)+2,ya,yb);sol=bvp4c(odefun,bcfun,sinit)sol=bvp4c(odefun,bcfun,sinit)t=linspace(0,4,101);t=linspace(0,4,101);y=y=deval(sol,tdeval(sol,t););plot(t,y(1,:),sol.x,sol.y(1,:),o,s
15、init.x,plot(t,y(1,:),sol.x,sol.y(1,:),o,sinit.x,sinit.y(1,:),s)sinit.y(1,:),s)legend(legend(解曲线解曲线,解点解点,粗略解粗略解) )刚性方程组刚性方程组例题：求解方程组例题：求解方程组解：先写解：先写M M函数函数eg6_6fun.meg6_6fun.m%M%M函数函数eg6_6fun.meg6_6fun.mfunction f=eg6_6fun(t,y)function f=eg6_6fun(t,y)f(1)=-0.01*y(1)-99.99*y(2);f(1)=-0.01*y(1)-99.99*y
16、(2);f(2)=-100*y(2);f(2)=-100*y(2);f=f(:);f=f(:);clear;tic;t,yclear;tic;t,y=ode45(eg6_6fun,0,400,2,1);toc,plot(t,y);=ode45(eg6_6fun,0,400,2,1);toc,plot(t,y);Elapsed_timeElapsed_time=4.5320=4.5320刚性方程组刚性方程组n这种问题的特点是这种问题的特点是y2y2下降很快，而下降很快，而y1y1下降太慢。下降太慢。u一方面，由于一方面，由于y2y2下降太快，为了保证数值稳定性，步长下降太快，为了保证数值稳定性，
17、步长h h需足需足够小；够小；u另一方面，由于另一方面，由于y1y1下降太慢，为了反映解的完整性，时间区间下降太慢，为了反映解的完整性，时间区间需足够长，这就造成计算量太大。需足够长，这就造成计算量太大。n这种方程组称病态方程组，即刚性方程组。这种方程组称病态方程组，即刚性方程组。n用用ode15sode15s求解求解clear;tic;t,yclear;tic;t,y=ode45(eg6_6fun,0,400,2,1);toc,plot(t,y);=ode45(eg6_6fun,0,400,2,1);toc,plot(t,y);Elapsed_timeElapsed_time=0.0470=
18、0.0470建模实验：产品销售量的增长建模实验：产品销售量的增长n例题：电饭锅一类的家庭主妇购买的商品，其实物广例题：电饭锅一类的家庭主妇购买的商品，其实物广告的效果非常明显。经调查发现，电饭锅的销售速度告的效果非常明显。经调查发现，电饭锅的销售速度与当时的销量成正比。试建立数学模型以预测销量。与当时的销量成正比。试建立数学模型以预测销量。n模型模型1 1：close;fplot(exp(0.9*x),0,10);hold on;close;fplot(exp(0.9*x),0,10);hold on;n模型模型2 2：考虑市场容量的限制，市场最大需求考虑市场容量的限制，市场最大需求10001000万台万台t,xt,x=ode45(inline(0.9*x*(1-=ode45(inline(0.9*x*(1-x/1000),t,x),0,10,1);plot(t,x)x/1000),t,x),0,10,1);plot(t,x)n模型模型3 3 考虑产品的更新换代等考虑产品的更新换代等作业题：作业题：AppoloAppolo卫星的运动轨迹卫星的运动轨迹已知已知AppoloAppolo卫星的运动轨迹满足下面方程卫星的运动轨迹满足下面方程其中，其中，初值条件为初值条件为试求解并绘制试求解并绘制AppoloAppolo卫星轨迹图卫星轨迹图
