MATLAB及其在理工课程中的应用-常微分方程

常微分方程,主要内容,学习MATLAB有关求解常微分方程（组）的指令；刚性方程组问题；实际问题的微分方程建模及求解方法；,解常微分方程主要MATLAB指令ode45四、五阶Runge-kutta法ode23二、三阶Runge-kutta法ode113多步Adams算法odeset解ode选项设置ode23t适度刚性问题梯形算法ode15s刚性方程组多步Gear法ode23s刚性方程组二阶Rosebrock法ode23tb刚性方程组低精度算法bvpinit边值问题的预估解bvp4c边值问题解法deval微分方程解的求值,预备知识：常微分方程,微分方程的概念常微分方程偏微分方程线性常微分方程初等积分法(分离变量法、积分因子法、常数变易法、降阶法等)常系数线性微分方程一阶变系数常微分方程：叠加原理求相应齐次微分方程解+一个特解高阶常系数微分方程：求相应齐次微分方程基本解+常数变易法求特解初值问题数值解注：高阶常微分方程初值问题需先转化为一阶常微分方程组,一、初值问题求解,常用调用格式：t,y=ode45(odefun,tspan,y0)参数说明：odefun表示f(t,y)的函数句柄或inline函数，t是标量，y是标量或向量；tspan如果是二维向量t0,tf，表示自变量初值t0和tf；如果是高维向量t0,t1,tn，则表示输出结点列向量；y0表示初值向量y0；t表示结点列向量(t0,t1,tn)T;y表示数值解矩阵，每一列对应y的一个分量；完整调用格式：t,y=ode45(odefun,tspan,y0，options,p1,p2,)options为计算参数设置（如精度要求等），默认表示；p1，p2，附加传递参数，odefun表示为f(t,y,p1,p2,),一、初值问题求解,ode45是最常用的求解微分方程的指令。采用变步长四、五阶龙格-库塔法，适合高精度问题；ode23与ode45类似，但精度低一些；ode113采用多步法，高低精度均适合；注：此三指令对于刚性方程组不宜采用。ode23t，ode23s,ode23tb,ode15s都是求解刚性方程组的指令。,一、初值问题求解,例题：解微分方程odefun=inline(y-2*t/y,t,y);t,y=ode45(odefun,0,4,1);t,yplot(t,y,o-)ans=01.00000.05021.04900.10051.09590.15071.1408.3.95022.98394.00003.0006ode45(odefun,0,4,1);t,y=ode45(odefun,0:1:4,1);t,y,一、初值问题求解,例题：解微分方程组解：将变量x，y合写为向量变量x，做M函数eg6_3fun.m%M函数eg6_3fun.mfunctionf=eg6_3fun（t，x）f(1)=-x(1)3-x(2);f(2)=x(1)-x(2)3;f=f(:);clear;t,x=ode45(eg6_3fun,030,1;0.5);subplot(1,2,1);plot(t,x(:,1),t,x(:,2),:);subplot(1,2,2);plot(x(:,1),x(:,2);,一、初值问题求解,例题：解微分方程组（竖直加热板的自然对流）已知：当=0时，解：首先引入辅助变量将方程组化为一阶方程组,一、初值问题求解,先写M函数eg6_4fun.m%M函数eg6_4fun.mfunctionf=eg6_4fun(t,y)f=y(2);y(3);-3*y(1)*y(3)+2*y(2)2-y(4);y(5);-2.1*y(1)*y(5);y0=0,0,0.68,1,-0.5;t,y=ode45(eg6_4fun,05,y0);plot(t,y(:,1),t,y(:,4),:),二、边值问题解法,常微分方程一阶方程组边值问题MATLAB标准形式：调用格式Sinit=bvpinit(tinit,yinit)由在粗略结点tinit的预估解yinit生成粗略解网格sinitSol=bvp4c(odefun,bcfun,sinit)odefun是微分方程组函数，bcfun为边值条件函数，sol.x为求解结点，sol.y是y(t)的数值解Sx=deval(sol,ti)计算由bvp4c得到的解在ti的值,二、边值问题解法,例题：求解边值问题解：首先改写成方程组边界条件为,二、边值问题解法,求解用M函数eg6_5fun.m%M函数eg6_5fun.mclear;close;sinit=bvpinit(0:4,1;0)odefun=inline(y(2);-abs(y(1),t,y);bcfun=inline(ya(1);yb(1)+2,ya,yb);sol=bvp4c(odefun,bcfun,sinit)t=linspace(0,4,101);y=deval(sol,t);plot(t,y(1,:),sol.x,sol.y(1,:),o,sinit.x,sinit.y(1,:),s)legend(解曲线,解点,粗略解),刚性方程组,例题：求解方程组解：先写M函数eg6_6fun.m%M函数eg6_6fun.mfunctionf=eg6_6fun(t,y)f(1)=-0.01*y(1)-99.99*y(2);f(2)=-100*y(2);f=f(:);clear;tic;t,y=ode45(eg6_6fun,0,400,2,1);toc,plot(t,y);Elapsed_time=4.5320,刚性方程组,这种问题的特点是y2下降很快，而y1下降太慢。一方面，由于y2下降太快，为了保证数值稳定性，步长h需足够小；另一方面，由于y1下降太慢，为了反映解的完整性，时间区间需足够长，这就造成计算量太大。这种方程组称病态方程组，即刚性方程组。用ode15s求解clear;tic;t,y=ode45(eg6_6fun,0,400,2,1);toc,plot(t,y);Elapsed_time=0.0470,建模实验：产品销售量的增长,例题：电饭锅一类的家庭主妇购买的商品，其实物广告的效果非常明显。经调查发现，电饭锅的销售速度与当时的销量成正比。试建立数学模型以预测销量。模型1：close;fplot(exp(0.9*x),0,10);holdon;模型2：考虑市场容量的限制，市场最大需求10