MATLAB解常微分方程及其应用

MATLAB解常微分方程及其应用授课教师：唐跃龙Telephone:QQ:

794218136Email:1．微分方程的根本概念定义1含有未知函数的导数(或微分)，同时也可能含有未知函数与自变量的方程，称为微分方程，未知函数是一元函数的，称为常微分方程；未知函数是多元函数的，称为偏微分方程.定义2微分方程中所出现的未知函数的最高阶导(或微分)的阶数，叫做微分方程的阶.定义3假设是的一次有理整式，称为阶线性常微分方程.定义4不是线性的微分方程称为非线性常微分方程.定义8假设某函数代入微分方程能使该方程成为恒等式，该函数就叫做该微分方程的解.定义5不含有未知函数及其导数(或微分)的项称为微分方程自由项.定义6自由项为零的微分方程称为齐次方程.定义7自由项不为零的微分方程称为非齐次方程.例1.指出以下微分方程的阶数并判断是线性还是非线性、齐次还是非齐次方程：定义5确定了通解中任意常数以后得到的解称为微分方程的特解.定义4假设微分方程的解中所含相互独立的任意常数的个数与微分方程的阶数相同，这样的解叫做微分方程的通解.例2.验证函数

是微分方程

的解.例3.求方程

的通解.例4.求方程

满足初始条件：当

时

的特解.注：这里所求的微分方程的解、通解、特解都是解析解〔精确解〕.2．MATLAB求微分方程的解析解格式一：y=dsolve('f1','f2',…,'fm')格式二：y=dsolve('f1','f2',…,'fm'

,'x')%指明自变量

fi既可以描述微分方程，又可描述初始条件或边界条件。如：描述微分方程时描述条件时例5求方程的解析解.解：>>symsx,t;x=dsolve('D1x=x*(1-x^2)','t')x=10-1(-1/(exp(C3-2*t)-1))^(1/2)假设方程变为：解：>>symsx,t;x=dsolve('Dx=x*(1-x^2)+1')警告:Explicitsolutioncouldnotbefound;implicitsolutionreturned.因此，不是所有方程都有解析解，只有极少局部的微分方程有解析解.当找不到原方程的解析解时，我们只能考虑其数值解.以下面的一阶常微分方程(ODE)初值问题为例：

数值解法就是求y(x)在某些分立的节点xn

上的近似值yn，用以近似y(xn)3.1简单欧拉(L.Euler,1707-1783)方法。

欧拉数值算法就是由初值通过递推求解，递推求解就是从初值开始，后一个函数值由前一个函数值得到。关键是构造递推公式。3．欧拉近似方法欧拉数值算法递推公式构造3.1.1差分法差分法就是用差商近似代替微商，即代入微分方程得到：对于等间隔节点可以得到：在xn节点上，微分方程可以写为作如下近似：那么得到欧拉解法递推公式的一般形式：具体求解过程为：简单欧拉方法程序function[outx,outy]=MyEuler(fun,x0,xt,y0,PointNum)%MyEuler用前向差分的欧拉方法解微分方程%fun表示f〔x，y〕%x0,xt表示自变量的初值和终值%y0表示函数在x0处的值,其可以为向量形式%PointNum表示自变量在[x0,xt]上取的点数ifnargin<5|PointNum<=0%如果函数仅输入4个参数值，那么PointNum默认值为100PointNum=100;endifnargin<4%y0默认值为0y0=0;endh=(xt-x0)/PointNum;%计算步长hx=x0+[0:PointNum]'*h;%自变量数组y(1,:)=y0(:)';%将输入存为行向量，输入为列向量形式fork=1:PointNumf=feval(fun,x(k),y(k,:));%计算f(x,y)在每个迭代点的值f=f(:)';y(k+1,:)=y(k,:)+h*f;%对于所取的点x迭代计算y值endouty=y;outx=x;%plot(x,y)%画出方程解的函数图例题：3.1.2折线法〔几何意义〕从〔x0,y0〕作曲线y(x)的切线，得切线方程：此切线与x=x1交点纵坐标为：y1yy(t)x0

x1x2

x3

x4

x5x6x从〔x1,y1〕作曲线y(x)在(x1,y(x1))的切线的平行线，得直线方程：与x=x2交点纵坐标为：y1yy(t)y2x1x2x3x4x5x6x7按照相似的方法，从〔xn,yn〕作曲线y(x)在(xn,y(xn))的切线的平行线，得直线方程：与x=xn+1交点纵坐标为：Q1折线近似曲线，n增大，误差变大x1x2x3x4x5x6x7yy(t)3.1.3数值积分对微分方程作从x0到x积分得：在[x0,x1]，积分采用矩形近似，得：同样，在[x0,x2]，积分采用矩形近似，得：同样，在[x0,xn+1]，积分采用矩形近似，得：作如下近似

得：3.1.4欧拉法误差利用泰勒级数得：作如下近似局部截断误差由泰勒展式可以看出，在欧拉方法中，假定yn

是精确的，由yn计算yn+1所产生的误差的数量级为O(h2)

。整体截断误差如果再考虑到yn本身的误差，计及误差的累积效果和局部截断误差，这种误差称为欧拉方法整体截断误差，数量级为O(h)

。

欧拉法是数值方法解微分方程最简单的一种，精度不高，实际应用不多，但它可以反映数值求解微分方程的特征。2.2改进的欧拉近似方法节点：步长：对于一般的常微分方程：积分法求欧拉公式时，矩形法采用前一点函数值求积，假设利用后一点，即为向后的欧拉方法。向后的欧拉方法递推公式为即向后的欧拉方法(隐式方法):预报---校正法

用欧拉方法预报用向后的欧拉方法校正

用欧拉方法预报用向后的欧拉方法校正---迭代积分法求欧拉公式时，积分采用梯形近似，即得到改进的欧拉方法。3.2.2改进的欧拉方法改进的欧拉方法递推公式为隐式方法:预报---校正法用欧拉方法预报用改进的欧拉方法校正计算公式为：例题：4．几种方法的比较1欧拉方法对于一般的常微分方程：取前一点的斜率作为平均斜率，即：整体截断误差：O(h)

2向后的欧拉方法取后一点的斜率作为平均斜率，即：整体截断误差：O(h)

3改进的欧拉方法取两点斜率平均值，即：整体截断误差：O(h2)

欧拉方法〔前、后〕和改进的欧拉方法优点是算法简单，但计算精度较低，不能满足实际问题求解对精度的要求，所以使用较少。

从以上三种方法的比较，可以推测：假设取多点处斜率的加权平均作为平均斜率，精度会更高，这就是龙格—库塔法的根本思想。四阶龙格—库塔法计算公式为：整体截断误差：O(h4)

例题：time_Euler=0.2500time_EulerPro=0.5000time_RK4=1.01505．MATLAB的常微分方程函数ode45ode23ode113ode15sode23sode23tode23tb格式[x,y]=ode45(′fun′,[x0,xn],y0,option]说明：适用于求解一阶常微分方程组fun定义微分方程组的函数文件名[x0,xn]求解区域y0初始条件向量option可选参数，由ODESET函数设置，比较复杂x输出自变量向量，y输出[y,y′,y″,..]没有一种算法可以有效地解决所有的ODE问题，因此MATLAB提供了多种ODE函数。例题二、三阶龙格库塔法例题四、五阶龙格库塔法t1=0.0310t2=0.0160例题：用MATLAB的符号解法,求解常微分方程：>>dsolve('Dy+3*x*y=x*exp(-x^2)')ans=(1/3*exp(-x*(x-3*t))+C1)*exp(-3*x*t)

>>dsolve('Dy+3*x*y=x*exp(-x^2)','x')ans=exp(-x^2)+exp(-3/2*x^2)*C1例题：用MATLAB的符号解法,求解常微分方程的特解：>>dsolve('x*Dy+2*y-exp(x)=0','y(1)=2*exp(1)','x')

ans=(exp(x)*x-exp(x)+2*exp(1))/x^2例题：采用ODE45求解描述某刚性问题的微分方程：functiondy=rigid(t,y)dy=zeros(3,1);%行向量dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);options=odeset('RelTol',1e-4,'AbsTol',[1e-41e-41e-5]);[T,Y]=ode45(@rigid,[012],[011],options);plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')legend('y1','y2','y3')例题：用MATLAB函数ode23,ode45,ode113,求解多阶常微分方程：functiondy=myfun03(x,y)dy=zeros(3,1)%初始化变量dydy(1)=y(2);%dy(1)表示y的一阶导数,其等于y的第二列值dy(2)=y(3);%dy(2)表示y的二阶导数dy(3)=2*y(3)/x^3+3*y(2)/x^3+3*exp(2*x)/x^3%dy(3)表示y的三阶导数%用ode23ode45ode113解多阶微分方程clear,clc[x23,y23]=ode23('myfun03',[1,10],[11030]);[x45,y45]=ode45('myfun03',[1,10],[11030]);[x113,y113]=ode113('myfun03',[1,10],[11030]);figure(1)%第一幅图plot(x23,y23(:,1),'*r',x45,y45(:,1),'ob',x113,y113(:,1),'+g')%作出各种函数所得结果legend('ode23解','ode45解','ode113解')title('ODE函数求解结果')figure(2)plot(x45,y45)%以ode45为例作出函数以及其各阶导数图legend('y','y一阶导数','y两阶导数')title('y,y一阶导数,y二阶导数函数图')例题：MATLAB在导热计算传热过程涉及面广，数学模型复杂。计算过程中涉及到许多运算方法。导热虽是容易做数学处理的一种热量传递方式，但其过程往往涉及常微分、偏微分方程、线性〔非线性〕方程组的求解。有一外径为4cm，内径为1.5cm，载有电流密度I为5000A/cm2的内冷钢质导体。导体单位时间发出的热量等于流体同时带走的热量，导体内壁面的温度为70℃。假定外壁面完全绝热。试确定，导体内部的温度分布；钢的导热系数k=0.38Kw/m·K，电导率ρ=2×1011kΩ·m。解：这是圆柱坐标中常物性一维稳态导热问题，结合此题具体条件导热微分方程式可简化为：结合边界条件，可得这一导热问题的数学描述为：此常微分方程的分析解，可调用MATLAB符号工具箱中的dsolve函数