数学实验常微分方程.doc

实验六 常微分方程的Matlab解法
一、实验目的
1． 了解常微分方程的解析解。
2． 了解常微分方程的数值解.
3． 学习掌握MATLAB软件有关的命令.
二、实验内容
一根长的无弹性细线，一段固定,另一端悬挂一个质量为的小球，在重力的作用下小球处于垂直的平衡位置.若使小球偏离平衡位置一个角度，让它自由,它就会沿圆弧摆动.在不考虑空气阻力的情况下,小球会做一定周期的简谐运动。利用牛顿第二定律得到如
下的微分方程

问该微分方程是线性的还是非线性的？是否存在解析解？如果不存在解析解,能否求出其近似解?
三、实验准备
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）
(2）
（3）
方程（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)
练习２ 求解微分方程

先求解析解，再求数值解，并进行比较。由
clear；
s=dsolve（'Dy=-y+t+1’，'y（0）=1’，'t’）
simplify(s）
可得解析解为.下面再求其数值解，先编写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 求方程

的数值解.不妨取。则上面方程可化为

先看看有没有解析解.运行MATLAB代码
clear;
s=dsolve（'D2y=9.8＊sin(y）’,'y(0）=15',’Dy(0)=0',’t')
simplify（s)
知原方程没有解析解。下面求数值解。令可将原方程化为如下方程组

建立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可见，随时间周期变化。
练习4 (刚性方程组求解）求下面刚性微分方程的解

使用dsolve可知解析解为

下面求数值解. 建立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给人的感觉似乎是始终大于0.5.但由的解析解可知,当时，两个分量均趋于0。下降极快,; 而下降很慢,（见下图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
可见计算太慢，需要48261步才能到达400.一方面,由于下降太快,为了保证数值稳定性，步长须足够小;另一方面,由于下降太慢,为了反映解的完整性，时间区间须足够长,这就造成计算量太大。这类方程称为刚性方程或病态方程。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吸引子） 求常微分方程

的数值解，初值取. 先建立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.

图6.7.5 Lorenz吸引子
实验作业
1．求下列微分方程的解析解
(1) 一阶线性方程
(2) 贝努利方程
(3) 高阶线性齐次方程
(4) 高阶线性非齐次方程
2．求方程

的解析解和数值解，并进行比较
3．分别用ode45和ode15s求解Van-del-Pol方程

的数值解,并进行比较.
4。 （Rosseler吸引子)用ode45数值求解方程

其中初值取,参数。
阅读资料
我国微分方程界的先辈申又枨教授
申又枨，数学家、数学教育家。从事函数论及微分方程的研究。主要成就涉及复变函数的插值理论。是在新中国建立微分方程学科研究的创始人之一。
　　申又枨，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年。