MATLAB解ODE

MATLAB 解ODE

分子二班 张雷 1314057

要求：

介绍MATLAB 求解常微分方程初值问题、刚性问题、边值问题，隐式微分方程、延迟微分方程、代数微分方程的功能函数和使用方法，并分别举例说明。含程序和运行结果。

常微分方程（ODE ）的定义：

凡含有参数，未知函数和未知函数导数 (或微分) 的方程，称为微分方程，有时简称为方程，未知函数是一元函数的微分方程称作常微分方程，未知数是多元函数的微分方程称作偏微分方程。微分方程中出现的未知函数最高阶导数的阶

数，称为微分方程的阶。定义式如下：

常微分方程的解析解：

有些微分方程可直接通过积分求解。有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解. 线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解。

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

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

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

常微分方程的数值解：

n )

,,",',()1()(-=n n y y y t f y )1(21,,',-===n n y y y y y y ?????????====-),,,,(''''211

3221n n n n y y y t f y y y y y y y

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

我们学过的有向前欧拉公式、向后欧拉公式、龙格-库塔公式、Adams公式等等等等。这其中，包含显格式、隐格式，单步法、多步法。

这些方法如果通过手工计算，费时费力，而且难免出错。使用Matlab进行计算，解决常微分方程问题就显示出优越性。

dsolve()函数

在MATLAB中，由函数dsolve()解决常微分方程（组）的求解问题，其具体格式如下：

r = dsolve('eq1,eq2,...', 'cond1,cond2,...', 'v')

'eq1,eq2,...'为微分方程或微分方程组，'cond1,cond2,...',是初始条件或边界条件，'v'是独立变量，默认的独立变量是't'。

函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。

求通解隐式微分方程

求解微分方程

输入语句：y=dsolve('Dy-x^2-y','x') (以字符串形式)

输出结果：y =C2*exp(x) - 2*x - x^2 – 2

边界问题

求解微分方程

输入语句：

syms x y;

y=dsolve('x*Dy+y-x=0','y(1)=2*exp(1)','x')

输出语句：

y =(x^2/2 + 2*exp(1) - 1/2)/x

求微分方程

输入语句：

syms x y;

y=dsolve('D2y-x=0','Dy(0)=2','y(0)=1','x')

输出语句：

y =x^3/6 + 2*x + 1

求解微分方程组

求

输入语句：[y1,y2]=dsolve('Dy1=4-2*y1','Dy2=2*x-y2','x')

输出语句；y1 =exp(-2*x)*(C10 + 2*exp(2*x))

y2 =exp(-x)*(C11 + 2*exp(x)*(x - 1))

问题四：初值问题

求解

输入语句：y=dsolve('Dy=4-2*x*sin(x)','y(0)=1','x')

输出语句：y =4*x - 2*sin(x) + 2*x*cos(x) + 1

数值解

ode23、ode45等功能函数

MATLAB提供了7个常微分方程求解器（solver），分别是ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb，其中前3个适用于求解非刚性（Nonstiff）问题，后4个适用于刚性问题。

基于龙格-库塔法，MATLAB提供了求常微分方程的数值解的函数。

一般调用格式为：

[t,y]=ode23(filename,tspan,y0)

[t,y]=ode45(filename,tspan,y0)

其中filename是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。,tspan形式为[t0,tf],表示求解区间，y0是初始状态列向量。t和y分别给出时间限量和相应的状态向量。

对于两者区别：ode是Matlab专门用于解微分方程的功能函数；solver有变步长（variable-step）和定步长（fixed-step）两种类型，不同类型有着不同的求解器。ode45求解器属于变步长的一种，采用Runge-Kutta算法；和他采用相同算法的变步长求解器还有ode23。ode45表示采用四阶，五阶runge-kutta单步算法,截断误差为(Δx)3。

初值问题

输入语句

function [ y ] = fun_01(x,y )

y=x^2+y

end

t0=0;tf=1; %确定t的范围

y0=2;

[t,y]=ode23('fun_01',[t0,tf],y0);

[T,Y] = ode45(odefun,tspan,y0,options)

options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等。

ode23、ode45都是变步长算法，但是也可以通过对tspan的设置，实现定步长计算。

还以上题为例，

输入:t=[0:0.5:1] ; %确定t的范围和步长y0=2;

[t,y]=ode23('fun_01',t,y0)

输出结果：

y = 2.0000

3.3448

5.8729

一些特殊情况

延迟微分方程：

延迟微分方程，数学领域中的一大模块。微分方程模型的共同特点是：系统状体变量（即未知函数）的导数不仅依赖于系统当前的状态，而且依赖于系统在过去某一时刻或某一历史时期的状态。

在数学领域中,时滞微分方程, 或延时微分方程(DDE) 是一类微分方程, 其中未知函数的在确定时刻的导数由先前时刻函数所决定.

微分代数方程：

微分代数方程就是几个微分方程和纯代数方程（没有导数）组成的一个系统。和偏微分方程类似，微分代数方程也很难找到精确的结果。

刚性问题：

常微分方程中同时包含有快变分量和慢变分量，二者变化速度相差非常大的量级，这种特点在数学上称为刚性。

描述这种过程的常微分方程组称为刚性方程组。

稳定步长h受绝对值最大的特征值 |λ|max控制，允许步长很小。

过程趋于稳定的时间由绝对值最小的特征值控制。但由于稳定性要求，仍要用小步长，计算量巨大，误差积累的影响也随着计算步数的增加越来越严重。

对于这些特殊情况，也可以通过Matlab进行求解

微分代数方程：

用矩阵形式表示出该DAEs

最后一个方程不是微分方程而是一个代数方程(这就是为什叫DAE的原因)，其实我们可以将它视为对三个状态变量的约束。

1.odefun=@(t,x)[-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);

2. 2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);

3. x(1)+x(2)+x(3)-1];%微分方程组

4. M=[1 0 0;0 1 0;0 0 0];%质量矩阵

5. options=odeset('mass',M);%对以DAE问题，mass属性必须设置

6. x0=[0.8;0.1;0.1];%初值

7. [t,x]=ode15s(odefun,[0 20],x0,options);%这里好像不能使用ode45

8. figure('numbertitle','off','name','DAE demo—by Matlabsky')

9. plot(t,x)

10.legend('x1(t)','x2(t)','x3(t)')

(其实，我真没看懂，只是复制了一下)

刚性问题：

MATLAB提供了7个常微分方程求解器（solver），分别是ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb，其中前3个适用于求解非刚性（Nonstiff）问题，后4个适用于刚性问题。

1.odefun=@(t,x)[0.04*(1-x(1))-(1-x(2))*x(1)+0.0001*(1-x(2))^2

2. -1e4*x(1)+3000*(1-x(2))^2];

3.x0=[0 1];

4.tspan=[0 100];

5.options=odeset('reltol',1e-6,'abstol',1e-8);

6.tic;[t2,y2]=ode15s(odefun,tspan,x0,options);toc

7.disp(['ode15s计算的点数' num2str(length(t2))])