实验七 常微分方程【实验目的】1． 了解常微分方程的基本概念。

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

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

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

【实验内容】如右图所示，一根长l 的无弹性细线，一段固定，另一端悬挂一个 质量为m 的小球，在重力的作用下小球处于垂直的平衡位置。

若使小球 偏离平衡位置一个角度θ，让它自由，它就会沿圆弧摆动。

在不考虑空气 阻力的情况下，小球会做一定周期的简谐运动。

利用牛顿第二定律得到如 下的微分方程0)0(',)0(,sin "0===θθθθθmg ml问该微分方程是线性的还是非线性的？是否存在解析解？如果不存在解析解，能否求出其近似解？【实验准备】1．微分方程的概念未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。

如果未知函数是一元函数，称为常微分方程。

常微分方程的一般形式为0),,",',,()(=n y y y y t F如果未知函数是多元函数，成为偏微分方程。

联系一些未知函数的一组微分方程组称为微分方程组。

微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。

若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为)()(')()(1)1(1)(t b y t a y t a y t a y n n n n =++++--若上式中的系数n i t a i ,,2,1),( =均与t 无关，称之为常系数或定常、自治、时不变的。

2．常微分方程的解析解有些微分方程可直接通过积分求解.例如,一解常系数常微分方程1+=y dtdy可化为dt y dy=+1,两边积分可得通解为1-=t ce y .其中c 为任意常数.有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解(显式解).线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。

高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。

一阶场微分方程与高阶微分方程可以互化，已给一个n 阶方程，),,",',()1()(-=n n y y y t f y设)1(21,,',-===n n y y y y y y ，可将上式化为一阶方程组⎪⎪⎪⎩⎪⎪⎪⎨⎧====-),,,,(''''2113221n n nn y y y t f y yy y y y y反过来，在许多情况下，一阶微分方程组也可化为高阶方程。

所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。

3．微分方程的数值解法除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解法。

考虑一阶常微分方程初值问题⎩⎨⎧=<<=000)()),(,()('y t y t t t t y t f t y f其中)'.,,,(,)',,,(,)',,,(020*******m m m y y y y f f f f y y y y ===所谓数值解法，就是寻求)(t y 在一系列离散节点f n t t t t ≤<<< 10上的近似值n k y k ,,1,0, =称k k k t t h -=+1为步长，通常取为常量h 。

最简单的数值解法是Euler 法。

Euler 法的思路极其简单：在节点出用差商近似代替导数ht y t y t y k k k )()()('1-≈+这样导出计算公式（称为Euler 格式）,2,1,0),,(1=+=+k y t hf y y k k k k他能求解各种形式的微分方程。

Euler 法也称折线法。

Euler 方法只有一阶精度，改进方法有二阶Runge-Kutta 法、四阶Runge-Kutta 法、五阶Runge-Kutta-Felhberg 法和先行多步法等，这些方法可用于解高阶常微分方程（组）初值问题。

边值问题采用不同方法，如差分法、有限元法等。

数值算法的主要缺点是它缺乏物理理解。

4．解微分方程的MATLAB 命令MATLAB 中主要用dsolve 求符号解析解，ode45,ode23,ode15s 求数值解。

ode23与ode45类似，只是精度低一些。

ode12s用来求解刚性方程组，是用格式同ode45。

可以用help dsolve, help ode45查阅有关这些命令的详细信息.练习1求下列微分方程的解析解（1）by+='ay（2）1=yy-xy=y''=,0)0(',)0(2sin()（3）1ff=ggf+fgg-)0('=)0(',,1='=,'方程（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 =5/3*sin(x)-1/3*sin(2*x)ans =-1/3*sin(x)*(-5+2*cos(x))方程（3）求解的MATLAB代码为:clear;s=dsolve('Df=f+g','Dg=g-f','f(0)=1','g(0)=1's =f: [1x1 sym]g: [1x1 sym]clear;s=dsolve('Df=f+g','Dg=g-f','f(0)=1','g(0)=1')simplify(s.f)simplify(s.g)s =f: [1x1 sym]g: [1x1 sym]ans =exp(t)*(cos(t)+sin(t))ans =exp(t)*(-sin(t)+cos(t))练习２求解微分方程=yt-yy)0(,1,1'=++先求解析解，再求数值解，并进行比较。

(1)求解析解clear;s=dsolve('Dy=-y+t+1','y(0)=1','t')simplify(s)s =t+exp(-t)ans =t+exp(-t)（2）求数值解%M函数 fun8.mfunction 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')ty3练习3 求方程0)0(',)0(,sin "0===θθθθθmg ml的数值解.不妨取15)0(,8.9,1===θg l .则上面方程可化为0)0(',15)0(,sin 8.9"===θθθθ（1）先求解析解 clear;s=dsolve('D2y=9.8*sin(y)','y(0)=15','Dy(0)=0','t') simplify(s)??? Error using ==> dsolveError, (in dsolve/IC) The 'implicit' option is not available when giving Initial Conditions. 知原方程没有解析解. （2）求数值解下面求数值解.令',21θθ==y y 可将原方程化为如下方程组⎪⎩⎪⎨⎧====0)0(,15)0()sin(8.9''211221y y y y y y%M 函数 fun9.mfunction f=fun9(t,y)f=[y(2),9.8*sin(y(1))]'; %f 向量必须为一列向量 命令窗口输入 clear;close;[t,y]=ode45('fun9',[0,10],[15,0]);plot(t,y(:,1)); %画θ随时间变化图,y(:2)则表示'θ的值xlabel('t'),ylabel('y1')0123456789101515.51616.5ty 1练习4 (刚性方程组求解)求下面刚性微分方程的解⎪⎩⎪⎨⎧==-=--=1)0(,2)0(,100'99.9901.0'2122211y y y y y y y （1）解析解clear;s=dsolve('Dy1=-0.01*y1-99.99*y2','Dy2=-100*y2','y1(0)=2','y2(0)=1','t') s =y1: [1x1 sym] y2: [1x1 sym]simplify(s.y1) simplify(s.y2)ans =exp(-100*t)+exp(-1/100*t) ans =exp(-100*t) （2）求数值解 %M 函数 fun10.mfunction f=fun10(t,y)f=[-0.01*y(1)-99.99*y(2),-100*y(2)]'; 命令窗口输入 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')012345678910-0.50.511.52tyclear;close;[t,y]=ode45('fun10',[0,400],[2,1]); tstep=length(t) minh=min(diff(t))maxh=max(diff(t))tstep =48261minh =5.0238e-004maxh =0.0102clear;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 =92minh =3.5777e-004maxh =32.1282050100150200250300350400-0.50.511.52ty[练习与思考]1．求下列微分方程的解析解。