数学实验常微分方程

实验六 常微分方程的Matlab 解法

一、实验目的

1． 了解常微分方程的解析解。

2． 了解常微分方程的数值解。

3． 学习掌握MATLAB 软件有关的命令。

二、实验内容

一根长l 的无弹性细线，一段固定，另一端悬挂一个质量为m 的小球，在重力的作用下小球处于垂直的平衡位置。若使小球偏离平衡位置一个角度θ，让它自由，它就会沿圆弧摆动。在不考虑空气阻力的情况下，小球会做一定周期的简谐运动。利用牛顿第二定律得到如 下的微分方程

0)0(',)0(,sin "0===θθθθθmg ml

问该微分方程是线性的还是非线性的？是否存在解析解？如果不存在解析解，能否求出其近似解？

三、实验准备

MATLAB 中主要用dsolve 求符号解析解，ode45,ode23,ode15s 求数值解。

ode45是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。ode23与ode45类似，只是精度低一些。ode12s 用来求解刚性方程组，是用格式同ode45。可以用help dsolve, help ode45查阅有关这些命令的详细信息.

四、实验方法与步骤

练习1 求下列微分方程的解析解

（1）b ay y +='

（2）1)0(',0)0(,)2sin(''==-=y y y x y

（3）1)0(',1)0(',','==-=+=g f f g g g f f

方程（1）求解的MATLAB 代码为：

clear;

s=dsolve('Dy=a*y+b')

结果为

s =-b/a+exp(a*t)*C1

方程（2）求解的MATLAB 代码为：

clear;

s=dsolve('D2y=sin(2*x)-y','y(0)=0','Dy(0)=1','x')

simplify(s) %以最简形式显示s

结果为

s =(-1/6*cos(3*x)-1/2*cos(x))*sin(x)+(-1/2*sin(x)+1/6*sin(3*x))*cos(x)+5/3*sin(x)

ans =-2/3*sin(x)*cos(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)

练习２ 求解微分方程

,1)0(,1'=++-=y t y y

先求解析解，再求数值解，并进行比较。由

clear;

s=dsolve('Dy=-y+t+1','y(0)=1','t')

simplify(s)

可得解析解为t e t y -+=。下面再求其数值解，先编写M 文件fun8.m

%M 函数fun8.m

function f=fun8(t,y)

f=-y+t+1;

再用命令

clear; close; t=0:0.1:1;

y=t+exp(-t); plot(t,y); %化解析解的图形

hold on; %保留已经画好的图形，如果下面再画图，两个图形和并在一起

[t,y]=ode45('fun8',[0,1],1);

plot(t,y,'ro'); %画数值解图形，用红色小圈画

xlabel('t'),ylabel('y')

结果见图6.7.1

图6.7.1 解析解与数值解

由图7.1可见，解析解和数值解吻合得很好。

下面我们讨论实验引例中的单摆问题.

练习3 求方程

0)0(',)0(,sin "0===θθθθθmg ml

的数值解.不妨取15)0(,8.9,1===θg l .则上面方程可化为

0)0(',15)0(,sin 8.9"===θθθθ

先看看有没有解析解.运行MATLAB 代码

clear;

s=dsolve('D2y=9.8*sin(y)','y(0)=15','Dy(0)=0','t')

simplify(s)

知原方程没有解析解.下面求数值解.令',21θθ==y y 可将原方程化为如下方程组

⎪⎩⎪⎨⎧====0)0(,15)0()

sin(8.9''2

1122

1y y y y

y y

建立M 文件fun9.m 如下

%M 文件fun9.m

function f=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')

结果见图6.7.2

图6.7.2 数值解

由图6.7.2可见，θ随时间t 周期变化。

练习4 (刚性方程组求解)求下面刚性微分方程的解

⎪⎩⎪⎨⎧==-=--=1)0(,2)0(,

100'99.9901.0'21

22211y y y y y y y 使用dsolve 可知解析解为

)100ex p(),

100ex p()01.0ex p(21t y t t y -=-+-=

下面求数值解. 建立M 文件fun10.m 如下

%M 文件fun10.m

function f=fun10(t,y) f=[-0.01*y(1)-99.99*y(2), -100*y(2)]';

运行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')

结果见图6.7.3

图6.7.3 数值解

图6.7.3给人的感觉似乎是1y 始终大于0.5.但由21,y y 的解析解可知,当∞→t 时,两个分量21,y y 均趋于0.2y 下降极快,0001.0)1.0(2

6.7.4).若用

clear; close;

[t,y]=ode45('fun10',[0,400],[2,1]);

tstep=length(t) %求计算总步数

minh=min(diff(t)) %最小步长

maxh=max(diff(t)) %最大步长

结果为

tstep =48261