MATLAB实验4-微分方程求解分析
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验4
实验名称：微分方程求解 实验目的：熟悉并掌握常微分方程数值求解的基 本原理、方法；掌握常微分方程模型的MATLAB 解析解法与数值解法。
实验任务： 1、P139，习题6 2、P140，习题13.
附录一：基于MATLAB的微分方程求解
一、微分方程的解析运算
1、函数求导的解析运算
y=diff(fun,x)
[x,y]=dsolve('D2x+2*Dx=x+2*y-exp(t)','Dy=4*x+3*y+4*exp(-t)')
二、微分方程数值求解
1、Euler法，梯形公式，预估-校正法，线性多 步法，Runge-Kutta等。 2、数值求解的MATLAB函数 (1) 、 ode45() ，是非刚性问题大部分场合的首选算 法。是单步算法，采用 4 ， 5 阶 Runge-Kutta 方程， 也称为RKF格式。 (2) 、 ode23() ，适用于非刚性问题精度较低的情 形，是单步算法，采用 2 ， 3 阶 Runge-Kutta 方程。
plot3(x(:,1),x(:,2),x(:,3));title(‘相空间三维图');
axis([10 42 -20 20 -20 25])
三、特殊微分方程的数值求解
1、刚性微分方程的解法 振动、电路及化学反应中出现刚性现象的线性常 系数方程组可表为
(t ) Ax f (t ) x
例：设著名的Lorenz模型的状态方程表示为
1 (t ) x1 (t ) x2 (t ) x3 (t ) x 2 (t ) x2 (t ) x3 (t ) x x 3 (t ) x1 (t ) x3 (t ) x2 (t ) x3 (t )
其来自百度文库，设
8 / 3, 10, 28.
其初始值为 ，求解该微分方程组。
x1 (0) x2 (0) 0, x3 (0) 1010
function xdot=lorenzeq(t,x)
xdot=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);x(1)*x(2)+28*x(2)-x(3)]; clear;x0=[0;0;1e-10]; [t,x]=ode45(@lorenzeq,[0,100],x0); subplot(1,2,1) plot(t,x); title(‘状态变量的时间响应图'); subplot(1,2,2)
syms t y u=exp(-5*t)*cos(2*t+1)+5;
uu=diff(u,t,2)+4*diff(u,t)+2*u;
dsolve('D4y+10*D3y+35*D2y+50*Dy+24*y=uu')
若考虑初始条件： y(0) 3, y(0) 2, y(0) y(3) (0) 0
y=diff(fun,x,n)
%函数fun对x求导数
%函数fun对x求n阶导数
例： f ( x) sin x /( x2 4x 3) 求 d 4 f ( x) / dx4 syms x f=sin(x)/(x^2+4*x+3);f1=diff(f,x,4);pretty(f1)
2、微分方程的解析运算
(3)、ode113()，适用于非刚性问题，采用Adams多步法。
(4)、ode23t()，适用于适度刚性问题，采用梯形算法。
(5)、ode15s()，适用于刚性问题，若ode45失败时可 尝试使用，采用多步法，Gear’s反向数值微分，精度 中等。
(6)、ode23s()，适用于刚性问题，采用单步法，2阶 Rosebrock算法，精度低。
常用的控制参数为： RelTol~相对误差容许上限，默认值为0.001， 为了保证较高的精度可以减少该值 。 AbsTol~为一个向量，其分量表示每个状态变量 允许的绝对误差，其默认值为10^(-6)，可以改变 求解精度。 例如，相对误差要求为10^(-7)时
options=odeset(‘RelTol’,1e-7);
下面，以ode45为例说明函数的用法。调用格式：
[t,x]=ode45(Fun,ts,x0)
[t,x]=ode45(Fun,ts,x0,options) 其中，Fun是由微分方程或方程组写成的M文件，ts描述方 程的求解区间，如果是二维向量[t0,tf]表示自变量初值t0和 终值 tf ，如果是高维向量 [t0,t1,…,tf]( 或者 t0:k:tf), 则表示输 出这些节点的值。x0为函数的初值，options为控制选项， 可以通过odeset()函数获取。
MATLAB 符号运算工具箱提供了一个线性常系数微分方 程的求解函数dsolve()，该函数允许用字符串的形式描述 微分方程及初值、边值条件，最终得出微分方程的解析 解。调用格式： y=dsolve(f1,f2, …,fm,’t’) 其中fi既可以描述微分方程，又可以描述初始条件或边 界条件。
y ( 4) (t ) 用D4y 表示，已知条件
则：
y=dsolve('D4y+10*D3y+35*D2y+50*Dy+24*y=uu','y( 0)=3','Dy(0)=2','D2y(0)=0','D3y(0)=0')
例2：求下列线性微分方程组
(t ) 2 x (t ) x(t ) 2 y (t ) e t x t y ( t ) 4 x ( t ) 3 y ( t ) 4 e
其中，x，f是n维向量，A是n阶矩阵。 k 为A的特征 根。称
max | Re(k ) | s min | Re(k ) |
为刚性比，s>10的方程可认为是刚性方程。
D2y(2)=3表示。
(2) 3 y
,可用
例1： u e5t
( 4)
cos(2t 1) 5 ，求下列微分方程的通解
( 3)
(t ) 50y (t ) 24y(t ) y (t ) 10y (t ) 35 y (t ) 4u (t ) 2u(t ) u
实验名称：微分方程求解 实验目的：熟悉并掌握常微分方程数值求解的基 本原理、方法；掌握常微分方程模型的MATLAB 解析解法与数值解法。
实验任务： 1、P139，习题6 2、P140，习题13.
附录一：基于MATLAB的微分方程求解
一、微分方程的解析运算
1、函数求导的解析运算
y=diff(fun,x)
[x,y]=dsolve('D2x+2*Dx=x+2*y-exp(t)','Dy=4*x+3*y+4*exp(-t)')
二、微分方程数值求解
1、Euler法，梯形公式，预估-校正法，线性多 步法，Runge-Kutta等。 2、数值求解的MATLAB函数 (1) 、 ode45() ，是非刚性问题大部分场合的首选算 法。是单步算法，采用 4 ， 5 阶 Runge-Kutta 方程， 也称为RKF格式。 (2) 、 ode23() ，适用于非刚性问题精度较低的情 形，是单步算法，采用 2 ， 3 阶 Runge-Kutta 方程。
plot3(x(:,1),x(:,2),x(:,3));title(‘相空间三维图');
axis([10 42 -20 20 -20 25])
三、特殊微分方程的数值求解
1、刚性微分方程的解法 振动、电路及化学反应中出现刚性现象的线性常 系数方程组可表为
(t ) Ax f (t ) x
例：设著名的Lorenz模型的状态方程表示为
1 (t ) x1 (t ) x2 (t ) x3 (t ) x 2 (t ) x2 (t ) x3 (t ) x x 3 (t ) x1 (t ) x3 (t ) x2 (t ) x3 (t )
其来自百度文库，设
8 / 3, 10, 28.
其初始值为 ，求解该微分方程组。
x1 (0) x2 (0) 0, x3 (0) 1010
function xdot=lorenzeq(t,x)
xdot=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);x(1)*x(2)+28*x(2)-x(3)]; clear;x0=[0;0;1e-10]; [t,x]=ode45(@lorenzeq,[0,100],x0); subplot(1,2,1) plot(t,x); title(‘状态变量的时间响应图'); subplot(1,2,2)
syms t y u=exp(-5*t)*cos(2*t+1)+5;
uu=diff(u,t,2)+4*diff(u,t)+2*u;
dsolve('D4y+10*D3y+35*D2y+50*Dy+24*y=uu')
若考虑初始条件： y(0) 3, y(0) 2, y(0) y(3) (0) 0
y=diff(fun,x,n)
%函数fun对x求导数
%函数fun对x求n阶导数
例： f ( x) sin x /( x2 4x 3) 求 d 4 f ( x) / dx4 syms x f=sin(x)/(x^2+4*x+3);f1=diff(f,x,4);pretty(f1)
2、微分方程的解析运算
(3)、ode113()，适用于非刚性问题，采用Adams多步法。
(4)、ode23t()，适用于适度刚性问题，采用梯形算法。
(5)、ode15s()，适用于刚性问题，若ode45失败时可 尝试使用，采用多步法，Gear’s反向数值微分，精度 中等。
(6)、ode23s()，适用于刚性问题，采用单步法，2阶 Rosebrock算法，精度低。
常用的控制参数为： RelTol~相对误差容许上限，默认值为0.001， 为了保证较高的精度可以减少该值 。 AbsTol~为一个向量，其分量表示每个状态变量 允许的绝对误差，其默认值为10^(-6)，可以改变 求解精度。 例如，相对误差要求为10^(-7)时
options=odeset(‘RelTol’,1e-7);
下面，以ode45为例说明函数的用法。调用格式：
[t,x]=ode45(Fun,ts,x0)
[t,x]=ode45(Fun,ts,x0,options) 其中，Fun是由微分方程或方程组写成的M文件，ts描述方 程的求解区间，如果是二维向量[t0,tf]表示自变量初值t0和 终值 tf ，如果是高维向量 [t0,t1,…,tf]( 或者 t0:k:tf), 则表示输 出这些节点的值。x0为函数的初值，options为控制选项， 可以通过odeset()函数获取。
MATLAB 符号运算工具箱提供了一个线性常系数微分方 程的求解函数dsolve()，该函数允许用字符串的形式描述 微分方程及初值、边值条件，最终得出微分方程的解析 解。调用格式： y=dsolve(f1,f2, …,fm,’t’) 其中fi既可以描述微分方程，又可以描述初始条件或边 界条件。
y ( 4) (t ) 用D4y 表示，已知条件
则：
y=dsolve('D4y+10*D3y+35*D2y+50*Dy+24*y=uu','y( 0)=3','Dy(0)=2','D2y(0)=0','D3y(0)=0')
例2：求下列线性微分方程组
(t ) 2 x (t ) x(t ) 2 y (t ) e t x t y ( t ) 4 x ( t ) 3 y ( t ) 4 e
其中，x，f是n维向量，A是n阶矩阵。 k 为A的特征 根。称
max | Re(k ) | s min | Re(k ) |
为刚性比，s>10的方程可认为是刚性方程。
D2y(2)=3表示。
(2) 3 y
,可用
例1： u e5t
( 4)
cos(2t 1) 5 ，求下列微分方程的通解
( 3)
(t ) 50y (t ) 24y(t ) y (t ) 10y (t ) 35 y (t ) 4u (t ) 2u(t ) u