内容发布更新时间 : 2024/6/10 18:01:54星期一 下面是文章的全部内容请认真阅读。
实验六 常微分方程的Matlab解法
一、实验目的
1． 了解常微分方程的解析解。 2． 了解常微分方程的数值解。
3． 学习掌握MATLAB软件有关的命令。
二、实验内容
一根长l的无弹性细线，一段固定，另一端悬挂一个质量为m的小球，在重力的作用下小球处于垂直的平衡位置。若使小球偏离平衡位置一个角度?，让它自由，它就会沿圆弧摆动。在不考虑空气阻力的情况下，小球会做一定周期的简谐运动。利用牛顿第二定律得到如 下的微分方程
ml?\?mgsin?,?(0)??0,?'(0)?0
问该微分方程是线性的还是非线性的？是否存在解析解？如果不存在解析解，能否求出其近似解？
三、实验准备
MATLAB中主要用dsolve求符号解析解，ode45,ode23,ode15s求数值解。
s=dsolve(‘方程1’, ‘方程2’,…,’初始条件1’，’初始条件2’ …,’自变量’) 用字符串方程表示，自变量缺省值为t。导数用D表示，2阶导数用D2表示，以此类推。S返回解析解。在方程组情形，s为一个符号结构。 [tout,yout]=ode45(‘yprime’,[t0,tf],y0) 采用变步长四阶Runge-Kutta法和五阶Runge-Kutta-Felhberg法求数值解，yprime是用以表示f(t,y)的M文件名，t0表示自变量的初始值，tf表示自变量的终值，y0表示初始向量值。输出向量tout表示节点(t0,t1, …,tn)T,输出矩阵yout表示数值解，每一列对应y的一个分量。若无输出参数，则自动作出图形。 ode45是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。ode23与ode45类似，只是精度低一些。ode12s用来求解刚性方程组，是用格式同ode45。可以用help dsolve, help ode45查阅有关这些命令的详细信息.
四、实验方法与步骤
练习1 求下列微分方程的解析解 （1）y'?ay?b
（2）y''?sin(2x)?y,y(0)?0,y'(0)?1 （3）f'?f?g,g'?g?f,f'(0)?1,g'(0)?1 方程（1）求解的MATLAB代码为：
clear;
s=dsolve('Dy=a*y+b')
结果为
1
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) 练习２ 求解微分方程
y'??y?t?1,y(0)?1,
先求解析解，再求数值解，并进行比较。由
clear;
s=dsolve('Dy=-y+t+1','y(0)=1','t') simplify(s)
可得解析解为y?t?e?t。下面再求其数值解，先编写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
1.41.351.31.25y1.21.151.11.051
00.20.40.60.81t 图6.7.1 解析解与数值解
2
由图7.1可见，解析解和数值解吻合得很好。 下面我们讨论实验引例中的单摆问题.
练习3 求方程
ml?\?mgsin?,?(0)??0,?'(0)?0
的数值解.不妨取l?1,g?9.8,?(0)?15.则上面方程可化为
?\?9.8sin?,?(0)?15,?'(0)?0
先看看有没有解析解.运行MATLAB代码
clear;
s=dsolve('D2y=9.8*sin(y)','y(0)=15','Dy(0)=0','t') simplify(s)
知原方程没有解析解.下面求数值解.令y1??,y2??'可将原方程化为如下方程组
??y1'?y2?y?2'?9.8sin(y1)y(0)?15,y(0)?0
?12建立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
16.5161y15.515012345678910
t 图6.7.2 数值解
由图6.7.2可见，?随时间t周期变化。
3
练习4 (刚性方程组求解)求下面刚性微分方程的解
?y1'??0.01y1?99.99y2? ?y2'??100y2,?y(0)?2,y(0)?12?1使用dsolve可知解析解为
y1?exp(?0.01t)?exp(?100t),下面求数值解. 建立M文件fun10.m如下
%M文件fun10.m
function f=fun10(t,y)
y2?exp(?100t)
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
21.51y1y0.5y20-0.501234
5t678910 图6.7.3 数值解
图6.7.3给人的感觉似乎是y1始终大于0.5.但由y1,y2的解析解可知,当t??时,两个分量
y1,y2均趋于0.y2下降极快,y2(0.1)?0.0001; 而y1下降很慢,y2(400)?0.0183(见下图
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
4
minh =5.0238e-004 maxh =0.0102
可见计算太慢,t需要48261步才能到达400.一方面,由于y2下降太快,为了保证数值稳定性,步长h须足够小;另一方面,由于y1下降太慢,为了反映解的完整性,时间区间须足够长,这就造成计算量太大.这类方程称为刚性方程或病态方程.ode45不适用于病态方程,下面我们用ode15s求解.
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 = 92
minh =3.5777e-004 maxh =32.1282
可见只需92步,最大步长为32,速度快了约500倍.函数图形见图6.7.4.
21.51y0.5y1y20-0.5050100150200250300350400
t 图6.7.4 数值解
练习5 (Lorenz吸引子) 求常微分方程
??dxdt??10x?10y???dy?28x?y?xz ?dt?dz??dt??83z?xy的数值解,初值取x(0)?y(0)?z(0)?1. 先建立M文件Lorenzf.m如下
%M文件Lorenzf.m
function f=lorenzf(t,x)
5