数学实验-常微分方程

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

一、实验目的

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

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

二、实验内容

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

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

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

三、实验准备

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

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''21

1221y 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

minh =5.0238e-004 maxh =0.0102

可见计算太慢,t 需要48261步才能到达400.一方面,由于2y 下降太快,为了保证数值稳定性,步长h 须足够小;另一方面,由于1y 下降太慢,为了反映解的完整性,时间区间须足够长,这就造成计算量太大.这类方程称为刚性方程或病态方程.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.

图6.7.4 数值解

练习5 (Lorenz 吸引子) 求常微分方程

?????????+-=--=+-=xy z dt

dz xz y x dt

dy

y x dt dx

38281010 的数值解,初值取1)0()0()0(===z y x . 先建立M 文件Lorenzf.m 如下

%M 文件Lorenzf.m

function f=lorenzf(t,x)

sig=10;bet=8/3;rho=28;

f=[sig*(x(2)-x(1)),x(1).*(rho-x(3))-x(2),x(1).*x(2)-bet*x(3)]; f=f(:);

运行MATLAB 代码

clear; close;

[t,y]=ode45('Lorenzf',[0,100],[1,1,1]); plot3(y(:,1),y(:,2),y(:,3)); xlabel x;ylabel y;zlabel z; 运行结果如图6.7.5.

x

y

z

图6.7.5 Lorenz 吸引子

实验作业

1．求下列微分方程的解析解

(1) 一阶线性方程2'3

=-y x y (2) 贝努利方程0'2=--y xy y

(3) 高阶线性齐次方程02'3"'"=+--y y y y (4) 高阶线性非齐次方程x y y y sin 32'3"=+- 2．求方程

3)0(',1)0(,'2")1(2===+y y xy y x

的解析解和数值解,并进行比较

3．分别用ode45和ode15s 求解Van-del-Pol 方程

()??

???===---1)0',0)0(0

)1(1000222x x x dt dx

x dt

x d 的数值解,并进行比较.

4. （Rosseler 吸引子）用ode45数值求解方程

?????????-+=+=+-=)()(c x z b dt

dz ay

x dt

dy

z y dt dx

其中初值取8)0(,3)0(,2)0(===z y x ，参数7.5,2.0,2.0===c b a .

阅读资料

我国微分方程界的先辈申又枨教授

申又枨，数学家、数学教育家。从事函数论及微分方程的研究。主要成就涉及复变函数的插值理论。是在新中国建立微分方程学科研究的创始人之一。

申又枨，1901年6月13日生于山西高平鼓楼。原名申祖佑，曾用名申幼声。后来由父亲改名为申又枨，其寓意是：在春秋战国时孔夫子的七十二个得意门生中有一位是山西人，名叫申又枨。父亲申声之，母亲李氏；夫人余嘉傲，生于1904年9月13日，曾任天津河北女子师范学院体育教师。长女申荔旋，次女申蔼旋，子申同健。

生平

1922至1926年在南开大学学习，一年级时是化学系学生，因对数学感兴趣，从二年级开始转到数学系学习；毕业后于1926至1927年期间在南开中学教书；申又枨于1927至1931年期间为南开大学助教；接着，于1931至1934年去哈佛大学数学系攻读博士学位，并在1935年得博士学位；1934至1935年在南开大学教书； 1935年应江泽涵教授的邀请，到北京大学数学系教书。抗日战争爆发后，随校共赴国难到昆明，是西南联大的教授。1945年抗战胜利后，随北京大学师生回到北平。在1947 至1949年期间，申又枨是北京大学数学系的代理系主任。1951 年他应邀去沈阳东北工学院数学系访问和工作。在1952年全国高等学校院系调整时，马寅初校长点名调申又枨教授回北大执教，并于1953年出任微分方程教研室首届主任，直到1978年4月22日逝世。

申又枨的一生经历了78个春秋，正逢中国历史多变的动荡时期，从封建的满清皇朝到军阀混战到中华民国，从抗日救亡运动到解放战争，从中华人民共和国的成立到社会主义的建设，从“文化大革命”的动乱到“四人帮”的覆灭，其中先生的许多喜怒哀乐犹如空中烟云俱往矣，惟倾心的事业及其献身精神将有传于世。

个人经历

1901年6月13日 出生于山西省高平鼓楼。 1922年8月-1926年7月 在天津南开大学学习。 1926年8月-1927年7月 在天津南开中学任教。

1927年8月-1931年7月在天津南开大学任教。

1931年8月-1934年7月在美国哈佛大学攻读数学博士。

1934年8月-1935年7月在天津南开大学任教。

1935年8月-1945年在西南联合大学任教。

1946年-1978年在北京大学任教。

主要成就

新中国微分方程学科的先驱

中华人民共和国成立之初，在党中央的号召下，当时知识分子全面学习苏联，把《理论联系实际》作为重要的学术标准，而且认真探索新路。当时数学界从苏联的经验中看到，微分方程和概率统计是数学联系实际的两大触角。这样，在全国高等学校院系调整后，北京大学先后于1953年和1955年设立了微分方程教研室和概率论教研室，并由系主任段学复邀请申又枨教授和许宝教授分别出任教研室主任。

其时，申又枨教授原来的研究方向是复变函数的插值理论，其中也涉及一些由微分方程定义的特殊函数。在解放前，他已经有意于对数学物理方程（或一些特殊的微分方程）的研究，曾打算研究空气动力学中超音速机翼的振动问题。这次被任命为微分方程教研室主任，他不辱使命，致力于教书育人的工作。

他在其职，谋其事，为教研室的发展而操劳。首先是帮助年轻教师提高教学水平，为此教研室从1953年开始同时举办了两个读书班：常微分方程方面攻读俄文版的斯捷潘诺夫的《微分方程教程》，而偏微分方程方面则是俄文版的吉洪诺夫的《数学物理方程》。这是大学本科二三年级的两本教科书，但当时对一些“提前毕业”的年轻助教而言，其中也有难点。例如，微分方程的解关于初值的可微性定理包含着不易理解的分析推导；又如，数理方程的波动理论涉及并不简单的物理概念，等等。申先生对我们的问题总是耐心听取，并且一起讨论，帮助解决。当不能立刻解答时，他就把问题带回家，在精心准备以后，再向我们作详细的讲解。这里不单是问题的解答，通常还带着学习方法的启示。

为祖国的数学发展默默奉献

1）申又枨教授在东北工学院教书期间（1950至1951学年）曾去哈尔滨访问。也许是地理位置的关系，哈尔滨是中华人民共和国成立后最早有俄文书店的城市。申又枨喜欢逛那里的书店，在一次偶然的浏览中发现了有两本使他爱不释手的俄文书。一本是B.涅梅茨基和B.斯捷潘诺夫的《微分方程定性理论》，另一本是索伯列夫的《泛函分析》，它们对当时的中国数学界还是陌生的著作。申又枨购买回沈阳后，就组织朋友把它们分别译成中文出版。后来，这两个中译本在国内成了流行的读物。事实上，申又枨教授非常重视数学新分支的萌芽。

2）1954年申又枨教授受国家教育部的委托，在北京大学举办了面向全国的暑期讲习班，为各地高等院校培训常微分方程和数学物理方程基础课的师资。常微分方程

和数理方程的课程大纲分别由申又枨教授和吴新谋教授主持制定。这在中国数学界属于最早的暑期讲习班之一。听讲者十分踊跃，挤满了北京大学第一教学楼的大教室，其中既有风华正茂的年轻人，也有白发苍苍的老教师，在炎热的夏天他们求知若渴。主讲者除资深教授申又枨、吴新谋和彭桓武外，还有年轻讲师谷超豪和叶彦谦等，他们也都怀着光荣的使命感，为国家培养急需的教学人才而不遗余力，其实当时的苏联教材刚引进新中国，他们也是边学边讲。教员与学员齐心合力，学习气氛非常浓厚，讲习班是名符其实，达到了预期的目的。

主要论著

一）Shen, Yu-cheng. Thesis, Harvard University 1936年 .

二）Shen, Yu-cheng. Interpolation to certain analytic functions by rational functions 1946年 .

三）Shen Yu-cheng. Interpolation to some calsses of analytic functions by rational functions with pre-assigned poles 1947年.