常微分方程课件

常微分方程

ORDINARYDIFFIERENTIAL

EQUATION(O.D.E)

1常微分方程【实验目的】

1.了解常微分方程的基本概念。2.了解常微分方程的解析解。3.了解常微分方程的数值解。4.学习、掌握MATLAB软件有关的命令。2

【实验内容】

如右图所示，一根长的无弹性细线，一端固定，另一端悬挂质量为的小球。在重力的作用下小球处于平衡位置。若使小球偏离平衡位置一定角度，放开它，它就会沿圆弧摆动。在不考虑空气阻力的情况下小球会做一定周期的简谐运动。利用牛顿第二定律得到如下的微分方程问该微分方程是线性的还是非线性的？是否存在解析解？如果不存在解析解，能否求出其近似解？3【实验准备】

1.微分方程的概念

未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为4如果未知函数是多元函数，称为偏微分方程。联系一些未知函数的一组微分方程称为微分方程组。微分方程中出现的未知函数的导数的最高阶数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般形式为若上式中的系数均与无关，称之为常系数或定常、自治、时不变的。52.常微分方程的解析解有些微分方程可直接通过积分求解。例如，一阶常系数微分方程可化为，两边积分可得通解为有些常微分方程求解可用一些技巧，如分离变量法，积分因子法，常数变异法，降阶法等化为可积分的方程求的解析解（显式解）。线性常微分方程的解满足叠加原理，从而他们的求解可归结为求解一个特解和相应齐次微分方程的通解。一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。

6一阶常微分方程组与高阶常微分方程可以互化，已知一个阶方程设，可将上式化为一阶常系数方程组反过来，在许多情况下，一阶常微分方程组也可化为高阶常微分方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。73.微分方程的数值解法除常系数线性微分方程可用特征根求解，少数特殊方程可用初等积分法求解外，大部分微分方程无解析解，应用中主要依靠数值解法。考虑一阶常微分方程初值问题其中所谓数值解法就是寻求在一系列离散节点上的近似值，称为步长，通常取为常量。最简单的数值解法是Euler法。

8Euler法的思想很简单：在节点处用差商近似代替导数这样导出计算公式（称为Euler格式）它能求解各种形式的微分方程。Euler法也称为折线法。

Euler法只有一阶精度，改进方法有二阶Runge-Kutta法、四阶Runge-Kutta法、五阶Runge-Kutta-Felhberg法和线性多步法等，这些方法可用于高阶常微分方程（组）初值问题。边值问题采用不同方法，如差分法、有限元法等。数值解法的主要缺点是它缺乏物理意义。94.解微分方程的MATLAB命令MATLAB中主要用dsolve求符号解析解，ode45,ode23,ode15s求数值解。S=dsolve(‘方程1’，‘方程2’，…’初始条件1’，‘初始条件2’，…‘自变量’)。用字符串表示方程，自变量缺省时为t.导数用D表示，2阶导数用D2表示，以此类推。S返回解析解。在方程组情形，S为一个符号结构。10[tout,yout]=ode45(‘yprime’,[t0,tf],y0)采用变步长四阶Rung-Kutta法和五阶Rung-Kutta-Felhberg法求数值解，yprime是表示f(t,y)的M文件名，t0表示自变量的初始值，tf表示自变量的终值，y0表示初始向量值。输出向量tout表示节点，输出矩阵yout表示数值解，每一列对应y的一个分量。若无输出参数，则自动作出图形。

ode45是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。ode23与ode45类似，只是精度低一些。ode15s用来求解刚性方程组，调用格式同ode45。可以用helpdsolve，helpode45查阅有关这些命令的详细信息。11【实验方法与步骤】练习1求下列微分方程的解析解（1）（2）（3）方程（1）求解的MATLAB代码为>>clear;>>s=dsolve('Dy=a*y+b')结果为s=-b/a+exp(a*t)*C112方程（2）求解的MATLAB代码为>>clear;>>s=dsolve('D2y=sin(2*x)-y','y(0)=0','Dy(0)=1','x')>>simplify(s)%以最简形式显示s结果为s=(-1/2*sin(x)+1/6*sin(3*x))*cos(x)+(-1/6*cos(3*x)-1/2*cos(x))*sin(x)+5/3*sin(x)ans=-2/3*cos(x)*sin(x)+5/3*sin(x)方程（3）求解的MATLAB代码为>>clear;>>s=dsolve('Df=f+g','Dg=g-f','f(0)=1','g(0)=1');>>simplify(s.f)%s是一个结构>>simplify(s.g)结果为ans=exp(t)*cos(t)+exp(t)*sin(t)ans=-exp(t)*sin(t)+exp(t)*cos(t)

13练习2求解微分方程先求解析解，再求数值解，并进行比较。由>>clear;>>s=dsolve('Dy=-y+t+1','y(0)=1','t');>>simplify(s)可得解析解为s=t+exp(-t)。下面再求其数值解，先编写M文件%M函数fun8.mfunctionf=fun8(t,y)y=-y+t+1;再运行相应的MATLAB代码：14clear;close;t=0:0.1:1;y=t+exp(-t);plot(t,y);%画解析解的图形holdon%保留已经画好的图形，如果下面要再画图，两个图形和并在一起[t,y]=ode45('fun8',[0,1],1);plot(t,y,'ro');%画数值解图形，用红色小圈画xlabel('t'),ylabel('y')运行结果见下图7.2。15图7.2解析解与数值解由图7.2可见，解析解和数值解吻合的很好。16下面我们讨论实验引例中的单摆问题练习3求方程的数值解。不妨设。则上面方程可化为先看看有没有解析解。MATLAB代码>>clear;>>s=dsolve('D2y=9.8*sin(y)','y(0)=15','Dy(0)=0','t')>>simplify(s)知原方程没有解析解。下面求数值解。令可将原方程化为如下方程组17建立M文件如下：%M文件fun9.mfunctionf=fun9(t,y)f=[y(2),9.8*sin(y(1))]‘;%f向量必须为一列向量MATLAB代码>>clear;>>close;>>[t,y]=ode45('fun9',[0,10],[15,0]);>>plot(t,y(:,1));%画随时间变化图，y(:,2)则表示的值>>xlabel('t'),ylabel('y1')运行结果见图318图7.3数值解由图7.3可见，随时间周期变化。图7.3数值解由图7.3可见，随时间周期变化。19练习4（刚性方程组求解）求下面刚性微分方程的解使用dsolve可知解析解为下面求数值解。建立M文件如下：%M文件fun10.mfunctionf=fun10(t,y)f=[-0.01*y(1)-99.99*y(2),-100*y(2)]'20运行MATLAB代码>>clear;close;>>[t,y]=ode45('fun10',[0,10],[2,1]);>>plot(t,y);text(1,1.1,'y1');text(1,0.1,'y2');>>xlabel('t');ylabel('y')运行结果见图7.4。

图7.4数值解21图7.4给人的感觉是始终大于0.5。但由的解析解可知，当时，两个分量均趋于0。下降极快，（0.1）<0.0001;而下降很慢，(400)0.0183(见下图7.5)。若用>>clear;close;>>[t,y]=ode45('fun10',[0,400],[2,1]);>>tstep=length(t)%求计算总步数>>minh=min(diff(t))%最小步长>>maxh=max(diff(t))%最大步长结果为tstep=48261minh=5.0238e-004maxh=0.0102可见计算太慢，需要48261步才能达到400。一方面，由于下降太快，为了保证数值稳定性，步长须足够小；另一方面，由于下降太慢，为了反映解的完整性，时间区间须足够长，这就造成计算量太大。22这类方程成为刚性方程或病态方程。Ode45不适用于病态方程,下面我们用Ode15s求解，相应的MATLAB代码为>>clear;close;>>[t,y]=ode15s('fun10',[0,400],[2,1]);>>plot(t,y);text(100,0.5,'y1');text(1,0.1,'y2');>>xlabel('t'),ylabel('y')>>tstep=length(t)>>minh=min(diff(t)>>maxh=max(diff(t)结果为tstep=92minh=3.5777e-004maxh=32.1282可见只需92步，最大步长为32，速度快乐约500倍。函数见图7.5。23图7.5数值解24【练习与思考】1求下列微分方程的解析解。a)一阶线性方程;b)贝努利方程;c)高阶线性齐次方程