实验五 用matlab求解常微分方程

1．微分方程的概念

未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为

F(t,y,y',y\,?,y(n))?0

如果未知函数是多元函数，成为偏微分方程。联系一些未知函数的一组微分方程组称为微分方程组。微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为

y(n)?a1(t)y(n?1)???an?1(t)y'?an(t)y?b(t)

若上式中的系数

ai(t),i?1,2,?,n均与t无关，称之为常系数。

dy 2．常微分方程的解析解

有些微分方程可直接通过积分求解.例如,一解常系数常微分方程dtdyy?1?dty?ce?1t?y?1可化为

,两边积分可得通解为.其中c为任意常数.有些常微分方程可用一些

技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解.

线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。

一阶常微分方程与高阶微分方程可以互化，已给一个n阶方程

y(n)?f(t,y',y\,?,y(n?1))

设

，可将上式化为一阶方程组 ?y1'?y2?y'?y3?2????y'?yn?n?1??yn'?f(t,y1,y2,?,yn)

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

3．微分方程的数值解法

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

?y'(t)?f(t,y(t)),t0?t?tf??y(t0)?y0

y1?y,y2?y',?,yn?y(n?1)其中

y?(y1,y2,?,ym)',f?(f1,f2,?,fm)',y0?(y10,y20,?,ym0)'.所谓数值解法，就

t?t1???tn?tfy,k?0,1,?,n是寻求y(t)在一系列离散节点0上的近似值k称

hk?tk?1?tk为步长，通常取为常量h。最简单的数值解法是Euler法。

Euler法的思路极其简单：在节点出用差商近似代替导数

y(tk?1)?y(tk)hy'(tk)?

这样导出计算公式(称为Euler格式)

yk?1?yk?hf(tk,yk),k?0,1,2,?

他能求解各种形式的微分方程。Euler法也称折线法。

Euler方法只有一阶精度，改进方法有二阶Runge-Kutta法、四阶Runge-Kutta法、五阶Runge-Kutta-Felhberg法和先行多步法等，这些方法可用于解高阶常微分方程(组)初值问题。边值问题采用不同方法，如差分法、有限元法等。数值算法的主要缺点是它缺乏物理理解。

4．解微分方程的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)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')

结果为

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

结果见图7.1

1.41.351.31.251.21.151.11.051y00.20.4t0.60.81 图16.6.1 解析解与数值解

由图16.6.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??y2'?9.8sin(y1)?y(0)?15,y(0)?02?1

建立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')

结果见图7.2

16.516y115.51501234 5t678910 图7.2 数值解

由图7.2可见，?随时间t周期变化。

习题

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

2．求方程

(1?x)y\?2xy',y(0)?1,y'(0)?32

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

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

的数值解,并进行比较.

?d2xdx2?x?0?2?1000(1?x)dt?dt?x(0)?0,x'?0)?1??