实验七__用matlab求解常微分方程.doc

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

一、实验目的：
1、熟悉常微分方程的求解方法，了解状态方程的概念；
2、能熟练使用dsolve函数求常微分方程（组）的解析解；
3、能熟练应用ode45\ode15s函数分别求常微分方程的非刚性、刚性的数值解；
4、掌握绘制相图的方法
二、预备知识：
1．微分方程的概念
未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为

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

若上式中的系数均与无关，称之为常系数。
2．常微分方程的解析解
有些微分方程可直接通过积分求解.例如,一解常系数常微分方程可化为,两边积分可得通解为.其中为任意常数.有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解.
线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。
一阶常微分方程与高阶微分方程可以互化，已给一个阶方程

设，可将上式化为一阶方程组

反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。
3．微分方程的数值解法
除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题

其中所谓数值解法，就是寻求在一系列离散节点上的近似值称为步长，通常取为常量。最简单的数值解法是Euler法。
Euler法的思路极其简单：在节点出用差商近似代替导数

这样导出计算公式（称为Euler格式）

他能求解各种形式的微分方程。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）
（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')
结果见图7.1

图16.6.1 解析解与数值解
由图16.6.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')
结果见图7.2

图7.2 数值解
由图7.2可见，随时间周期变化。
以上这些都是常微分方程的精确解法，也称为常微分方程的符号解。但是，我们知道，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB具有丰富的函数，我们将其统称为solver，其一般格式为：
[T,Y]=solver(odefun,tspan,y0)
该函数表示在区间tspan=[t0,tf]上，用初始条件y0求解显式常微分方程
solver为命令ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb之一，这些命令各有特点。我们列表说明如下：
求解器
ODE类型
特点
说明
ode45
非刚性
一步算法，4,5阶Runge-Kutta
方法累积截断误差
大部分场合的首选算法
ode23
非刚性
一步算法，2,3阶Runge-Kutta
方法累积截断误差
使用于精度较低的情形
ode113
非刚性
多步法，Adams算法，高低精度均可达到
计算时间比ode45短
ode23t
适度刚性
采用梯形算法
适度刚性情形
ode15s
刚性
多步法，Gear’s反向
数值积分，精度中等
若ode45失效时，
可尝试使用
ode23s
刚性
一步法，2阶Rosebrock算法，
低精度。
当精度较低时，
计算时间比ode15s短
odefun为显式常微分方程中的
tspan为求解区间，要获得问题在其他指定点上的解，则令（要求单调），
y0初始条件。
例5：求解常微分方程，，的MATLAB程序如下：fun=inline('-2*y+2*x*x+2*x')；[x,y]=ode23(fun,[0,0.5],1)
结果为：
x = 0,0.0400,0.0900,0.1400,0.1900,0.2400,0.2900,0.3400,0.3900,0.4400,0.4900,0.5000
y = 1.0000,0.9247,0.8434,0.7754,0.7199,0.6764,0.6440,0.6222,0.6105,0.6084,0.6154,0.6179
例6：求解常微分方程的解，并画出解的图形。
分析：这是一个二阶非线性方程，用现成的方法均不能求解，但我们可以通过下面的变换，将二阶方程化为一阶方程组，即可求解。
令：，，，则得到：

接着，编写vdp.m如下：
function fy=vdp(t,x)
fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)];
再编写m文件sy12_6.m如下：
y0=[1;0]
[t,x]=ode45(@vdp,[0,40],y0);
y=x(:,1);dy=x(:,2);
plot(t,y,t,dy)

三、实验内容：
1．求下列微分方程的解析解

2．求方程

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

的数值解,并进行比较.