数学实验课件 第9章9.2.pptx

9.2微分方程的数值解;除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无解析解，应用中主要依靠数值解法．考虑一阶常微分方程初值问题

其中所谓数值解法，就是寻求y(x)在一系列离散节点上的近似值yk，．称为步长，通常取为常量h．最简单的数值解法是Euler法．;9.2.1欧拉法;9.2.2龙格-库塔法;在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.m

functionf=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.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.m

functionf=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.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))