数学实验课件：常微分方程实验

常微分方程实验微分方程差不多是和微积分同时先后产生的，苏格兰数学家耐普尔创立对数的时候，就讨论过微分方程的近似解.

牛顿在建立微积分的同时，对简单的微分方程用级数来求解．瑞士数学家雅各布·贝努利、欧拉、法国数学家克雷洛、达朗贝尔、拉格朗日等人又不断地研究和丰富了微分方程的理论．常微分方程的形成与发展是和力学、天文学、物理学，以及其他科学技术的发展密切相关的.

数学的其他分支的新发展，如复变函数、李群、组合拓扑学等，都对常微分方程的发展产生了深刻的影响，当前计算机的发展更是为常微分方程的应用及理论研究提供了非常有力的工具．

牛顿研究天体力学和机械动力学的时候，利用了微分方程这个工具，从理论上得到了行星运动规律.法国天文学家勒维烈和英国天文学家亚当斯使用微分方程各自计算出那时尚未发现的海王星的位置.这些都使数学家更加深信微分方程在认识自然、改造自然方面的巨大力量．本章介绍利用MATLAB求解常微分方程的解析解和数值解.常微分方程的解析解9.1CONTENTS9.2微分方程的数值解9.2.1欧拉法9.2.2龙格-库塔法9.1常微分方程的解析解

定义9.1

含有未知函数的导数的方程称为微分方程.

如果未知函数是一元函数，称为常微分方程.

常微分方程的一般形式为.

由几个微分方程联立而成的方程组称为微分方程组．

定义9.2

微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶．若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为.

若上式中的系数

均与

x无关，称之为常系数.

二阶及二阶以上的微分方程称为高阶微分方程．

一阶常微分方程与高阶微分方程可以互化，已给一个n阶方程设

，

可将上式化为一阶方程组

反过来，在许多情况下，一阶微分方程组也可化为高阶方程.

所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解．

MATLAB中主要用dsolve求符号解析解，调用格式如下：

S=dsolve(eqn,cond)

求带初始条件的微分方程特解，若缺省cond，则求微分方程的通解.利用“diff”和“==”表示微分方程．例如diff(y,x)==y表示dy/dx=y．例9.1

求下列微分方程的解析解（1）（2）（3）解（1）>>symsy(t)ab

>>eqn=diff(y,t)==a*y+b;

>>dsolve(eqn)

ans=

-(b-C1*exp(a*t))/a可知微分方程的通解为

，其中

，C1是任意常数．（2）>>symsy(x)>>eqn=diff(y,x,2)==sin(2*x)-y;>>Dy=diff(y,x);>>cond=[y(0)==0,Dy(0)==1];>>dsolve(eqn,cond)ans=(5*sin(x))/3-sin(2*x)/3可知微分方程的特解为.（3）>>symsf(t)g(t)>>eqns=[diff(f,t)==f+g,diff(g,t)==g-f];>>cond=[f(0)==1,g(0)==1];>>S=dsolve(eqns,cond)S=包含以下字段的struct:g:[1×1sym]f:[1×1sym]>>S.fans=exp(t)*cos(t)+exp(t)*sin(t)>>S.gans=exp(t)*cos(t)-exp(t)*sin(t)可知微分方程组的特解为.9.2微分方程的数值解

除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无解析解，应用中主要依靠数值解法．考虑一阶常微分方程初值问题其中所谓数值解法，就是寻求y(x)在一系列离散节点

上的近似值yk

，

．称

为步长，通常取为常量h．最简单的数值解法是Euler法．9.2.1欧拉法

Euler法的思路极其简单：在节点处用差商近似代替导数

这样导出计算公式（称为Euler格式）

它能求解各种形式的微分方程.Euler法也称折线法．

Euler方法只有一阶精度，改进方法有二阶Runge-Kutta法、四阶Runge-Kutta法、五阶Runge-Kutta-Felhberg法和先行多步法等，这些方法可用于解高阶常微分方程（组）初值问题.边值问题采用不同方法，如差分法、有限元法等.

数值算法的主要缺点是它缺乏物理理解．9.2.2龙格-库塔法

欧拉公式易于计算，但精度不高，收敛速度慢．在实际应用中，我们采用龙格-库塔（Runge-Kutta）方法，基本思想为进一步可得在MATLAB中，利用ode23，ode45求微分方程数值解．

[t,y]=ode23(odefun,tspan,y0)

[t,y]=ode45(odefun,tspan,y0)

求微分方程组y′=f(t,y)从t0到tf的积分，初始条件为y0.其中tspan=[t0tf].解数组y中的每一行都与列向量t中返回的值相对应．

ode45是最常用的求解微分方程数值解的命令，采用四阶和五阶Runge-Kutta算法，是一种自适应步长（变步长）的常微分方程数值解法，对于刚性方程组不宜采用.ode23与ode45类似，采用二阶和三阶Runge-Kutta算法，只是精度低一些．ode45是解决数值解问题的首选方法．例9.2

求解微分方程先求解析解，再求数值解，并进行比较.解>>clear;>>symsy(t);>>eqn=diff(y,t)==-y+t+1;>>cond=y(0)==1;>>dsolve(eqn,cond)ans=t+exp(-t)可得解析解为

．下面再求其数值解，先编写M文件fun92.m．%M函数fun92.mfunctionf=fun92(t,y)f=-y+t+1;再用命令>>clear;t=0:0.1:1;>>y=t+exp(-t);plot(t,y);

%绘制解析解的图形>>holdon;

%保留已经画好的图形，如果下面再画图，两个图形和并在一起>>[t,y]=ode45(@fun92,[0,1],1);>>plot(t,y,'ro');

%绘制数值解图形>>xlabel('t'),ylabel('y')结果见图9-1．

由图9-1所示，解析解和数值解吻合得很好．图9-1解析解与数值解例9.3已知方程当

时，上面方程可化为求上面方程的解析解和数值解．解先求解析解，>>symsy(t)>>eqn=diff(y,t,2)==9.8*sin(t);>>Dy=diff(y,t);>>cond=[y(0)==15,Dy(0)==0];>>dsolve(eqn,cond)ans=(49*t)/5-(49*sin(t))/5+15可知方程的解析解为

.下面求数值解．令

可将原方程化为如下方程组建立函数文件fun93.m如下%M文件fun93.mfunctionf=fun93(t,y)f=[y(2),9.8*sin(y(1))]';

%f向量必须为一列向量运行MATLAB代码>>clear;close;>>[t,y]=ode45(@fun93,[0,10],[15,0]);>>plot(t,y(:,1));

%画

随时间变化图,y(:2)则表示

的值>>xlabel('t'),ylabel('y1')结果见图9-2．

由图9-2可见，

随时间t周期变化．图9-2数值解

例9.4

Lotka-Volterra方程，也即捕食者-猎物模型的一对一阶常微分方程

（9-1）

变量

x

和

y

分别计算猎物和捕食者的数量.当没有捕食者时，猎物数量将增加，当猎物匮乏时，捕食者数量将减少.使用初始条件x(0)=y(0)=20，使捕食者和猎物的数量相等．求当α=0.01和β=0.02时方程的解．解在MATLAB中，两个变量x和y可以表示为向量y中的两个值．同样，导数是向量yp中的两个值.当α=0.01和β=0.02时，方程组（9-1）可以表示为：

yp(1)=(1-alpha*y(2))*y(1)

yp(2)=(-1+beta*y(1))*y(2)MATLAB中的函数文件lotka.m来表示Lotka-Volterra方程，typelotkafunctionyp=lotka(t,y)%LOTKALotka-Volterrapredator-preymodel.%Copyright1984-2014TheMathWorks,Inc.yp=diag([1-.01*y(2),-1+.02*y(1)])*y;首先使用ode23在区间0<t<15中求解lotka中定义的微分方程．t0=0;tfinal=15;y0=[20;20];[t,y]=ode23(@lotka,[t0tfinal],y0);绘制两个种群数量对时间的图，见图9-3a．plot(t,y(:,1),'--',t,y(:,2))title('捕食者-猎物模型')xlabel('时间')ylabel('数量')legend('猎物','捕食者','Location','North')再使用

ode45

求解该方程组，得到相轨线图9-3b.[T,Y]=ode45(@lotka,[t0tfinal],y0);plot(y(:,1),y(:,2),'-',Y(:,1),Y(:,2),'--');title('相轨线')legend('ode23','ode45')

图9-3捕食者-猎物模型a)捕食者—猎物模型b）相轨线

图9.3a中实线表示的是微分方程的解，虚线表示的是解曲线的导数所描绘的曲线.可以看出，猎物数量减少时，捕食者的数量也在减少，捕食者的种群数量会随着猎物种群数量的增加而不断增加．猎物的种群数量增加时，捕食者数量也在增加，但是当捕食者达到一定的程度后，猎物又在不断减少.即猎物种群数量增加→捕食者种群数量增加→猎物种群数量减少→捕食者种群数量减少→猎物种群数量增加→再次循环．这种变化趋势反映了生态系统中普遍存在的负反馈调节．

图9.3b中，使用不同的数值方法求解微分方程会产生略微不同的答案.可以看出利用ode45函数得到的图形是平滑的，要比ode23函数得到的结果精度更高．

例9.5（罗伦兹吸引子与“蝴蝶效应”）吸引子在1963年由美国麻省理工学院的气象学