实验七 常微分方程

实验七 常微分方程

【实验目的】

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 dt

dy

可化为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 n

n y y y t f y y

y y y y y

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

3．微分方程的数值解法

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

⎩⎨

⎧=<<=0

00)()),(,()('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 法的思路极其简单：在节点出用差商近似代替导数

h

t 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）b

y+

='

ay

（2）1

=y

y

-

x

y

=

y

''=

,0

)0('

,

)0(

2

sin(

)

（3）1

f

f

=g

g

f

+

f

g

g

-

)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)