Matlab 常微分方程

作者：未知信息来源:未知 2006-1-26

字体大小:小 中 大网友评论条 进入论坛

m]function Yd=DYdt(t,Y)% t 一定是标量形式的自变量% Y 必须是列向量global G ME % 在函数中定义全局变量传递参数xy=Y(1:2)。-G*ME*xy/r^3]。（4）对微分方程进行解算global G ME % 在主程序中定义全局变量传递参数 1G=6。ME=5。...

5.14 常微分方程
5.14.2 ODE解算指令的使用演示
5.14.2.1 解算指令简洁格式使用示例

【 * 例 5.14.2.1-1 】采用最简洁格式的 ODE 文件和解算指令，研究围绕地球旋转的卫星轨道。

图 5.14.2.1-1 地球轨道卫星运动加速度

（1）问题的形成
轨道上运动的卫星，在 Newton 第二定律，和万有引力定律作用下，有。即，，而。里引力常数，是地球的质量。又假定卫星以初速度在处进入轨道。

（2）构成一阶微分方程组

令，则

( 5.14.2 .1-5)

初始条件为
( 5.14.2 .1-6)

（3）根据式 (5.14.2.1-5) 编写最简洁格式的 ODE 文件
[dYdt.m]
function Yd=DYdt(t,Y)
% t 一定是标量形式的自变量
% Y 必须是列向量
global G ME % 在函数中定义全局变量传递数
xy=Y(1:2);Vxy=Y(3:4); % 前两个元素是 " 位移量 " ，后两个是 " 速度量 " 。
r=sqrt(sum(xy.^2));
Yd=[Vxy;-G*ME*xy/r^3]; %Yd 必须按式 (5.14.2.1-5) 编写，是与 Y 同维的列向量。

（4）对微分方程进行解算
global G ME % 在主程序中定义全局变量传递参数 <1>
G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;
tspan=[t0,tf]; % 指定解算微分方程的时间区间
Y0=[x0;0;0;vy0]; % 按式 (5.14.2.1-6) 给定初值向量
[t,YY]=ode45('DYDt',tspan,Y0); %<8>
X=YY(:,1); % 输出 Y 的第一列是位移数据 x(t)
Y=YY(:,2); % 输出 Y 的第二列是位移数据 y(t)
plot(X,Y,'b','Linewidth',2); hold on
axis('image') % 保证 x 、 y 轴等刻度，且坐标框恰包容图形
[XE,YE,ZE] = sphere(10); % 产生单位球面数据
RE=0.64e7; % 地球半径
XE=RE*XE;YE=RE*YE;ZE=0*ZE; % 坐标纸上的地球平面数据
mesh(XE,YE,ZE),hold off % 绘地球示意图

图 5.14.2.1-2 卫星轨道

【 * 例 5.14.2.1-2 】上例中，程序间的参数（如 G 和 ME ）传送，是依靠全局变量形式实现的。一般说来，编写程序时，应尽量少用全局变量，以免引起混乱。本例演示参数如何在指令间直接传送。要实现参数直接传送，必须对上例中的程序进行修改，具体如下：

（1）重写 ODE 文件
[DYDt2.m]
function Yd=DYDt2(t,Y,flag,G,ME)
% flag 按 ODE 文件格式规定，必须是第三输入宗量。对它的赋值由 ode45 指令自动产生。
% 第 4 、 5 宗量是被传递的参数
switch flag
case '' % 按规定：这里必须是空串。在此为 " 真 " 时，完成以下导数计算。
X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V;-G*ME*X/r^3];
otherwise
error([ 'Unknown flag ''' flag '''.' ]);
end

（2）按以下方法修改上例第（4）步中的程序，并运行之，可得相同结果。
删去原主程序第 <1> 条指令；把原主程序的第 <8> 条指令改写为
[t,YY]=ode45('DYDt2',tspan,Y0,[],G,ME); % 第 4 宗量取缺省设置
% 第 5 、 6 宗量是被传递参数

5.14.2.2 解算指令较复杂格式的使用示例

【 * 例 5.14.2.2-1 】带事件设置的 ODE 文件及主程序编写演示。本例将以较高精度计算卫星经过近地点和远地点的时间，并在图上标志。

（1） ODE 文件的编写
[DYDt3.m]
function varargout=DYDt3(t,Y,flag,G,ME,tspan,Y0)
% DYDt3.m 供主程序调用的 ODE 函数文件
% 本文件自带三个子函数： f,fi,fev 。
% t, Y 分别是自变量和一阶函数向量，是最基本的输入宗量。
% flag 第三输入宗量，它专供解算指令（如 ode45 ）作调用通知。
% 在运行中，解算指令会根据需要向 flag 发不同的字符串。
% varargout 是 " 变 " 输出宗量。它由变维的元胞数组构成。每个元胞中可以存放指令
% 所产生的任意形式的数据。
switch flag
case '' % 必须用空串符。情况为 " 真 " 时，计算导数 dY/dt = f(t,Y) 。
varargout{1} = f(t,Y,G,ME); % 输出为一个元胞，容纳 f 子函数的一个输出 Yd
case 'init' % 必须用 'init' ，情况为 " 真 " 时，传送计算区间、初值、设置参数。
[varargout{1:3}] = fi(tspan,Y0); % 输出为三个元胞，容纳 fi 子函数的三个输出量
case 'events' % 必须用 'events' ，情况为 " 真 " 时，设置事件性质。
[varargout{1:3}] = fev(t,Y,Y0); % 输出三个元胞，容纳 fev 子函数的三个输出量
otherwise
error([ 'Unknown flag ''' flag '''.' ]);
end

% ------------------------------------------------------------------
function Yd = f(t,Y,G,ME)
% 计算导数子函数，被 " 父 " 函数 DYDt3 调用。
X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V; -G*ME*X/r^3];
% ------------------------------------------------------------------
function [ts,y0,options] = fi(tspan,Y0)
% 设置时间区间、初值、算法参数子函数，被 " 父 " 函数 DYDt3 调用。
ts=tspan;y0 = Y0;
options = odeset( 'Events' , 'on' , 'Reltol' ,1e-5, 'Abstol' ,1e-4);
% 开动 ode45 的 " 事件判断 " 功能，设置相对误差和绝对误差。

% ------------------------------------------------------------------
function [value,isterminal,direction] = fev(t,Y,Y0)
% 事件判断子函数，被 " 父 " 函数 DYDt3 调用。
dDSQdt = 2 * ((Y(1:2)-Y0(1:2))' * Y(3:4));
%dDSQdt 是离初始点的二乘距离 u 的时间导数 du/dt ，而 u=(x-x0)^2+(y-y0)^2 。
value = [dDSQdt; dDSQdt]; % 定义两个穿越 0 的事件
direction = [1; -1]; % 第一事件：以渐增方式穿越 0 。第二事件：以渐减方式穿越 0
isterminal = [1; 0]; % 第一事件发生后，终止计算；而第二事件发生后，继续计算。

（2）运行以下主程序
G=6.672e-11;ME=5.97e24;vy0=4000; x0=-4.2e7;t0=0;tf=60*60*24*9;
tspan=[t0,tf];Y0=[x0;0;0;vy0];
[t,YY,Te,Ye,Ie]=ode45('DYDt3',[],[],[],G,ME,tspan,Y0); % <3>
X=YY(:,1);Y=YY(:,2);
plot(X,Y,'b','Linewidth',2);hold on
text(0,6e7,' 轨道 ','Color','b') % 产生蓝色文字注释
axis('image'); % 保证 x 、 y 轴等长刻度，且坐标框恰包容图形

% 在三个事件发生点上画标记
plot(Ye(1,1),0.4e7+Ye(1,2),'r^','MarkerSize',10)
plot(Ye(2,1),0.4e7+Ye(2,2),'bv','MarkerSize',10)
plot(Ye(3,1),-0.4e7+Ye(3,2),'b^','MarkerSize',10)

% 把轨道的半周期和全周期标在图上
text(0.8*Ye(3,1),-2e7+Ye(3,2),['t3=' sprintf('%6.0f',Te(3))])
text(0.8*Ye(2,1),1.5e7+Ye(2,2),['t2=' sprintf('%6.0f',Te(2))])

% 在 x-y 坐标上画地球
[XE,YE,ZE] = sphere(10);RE=0.64e7;XE=RE*XE;YE=RE*YE;ZE=0*ZE;
mesh(XE,YE,ZE)
text(1e7,1e7,' 地球 ','Color','r'), hold off % 产生红色文字注释

图 5.14.2.2-1 带事件标注的卫星轨道图

5.14.3 关于ODE文件的说明

（1） ODE 文件的模板
下面就是 MATLAB 的模板文件 odefile.m ，为便于读者阅读，本书作者适当地给以注解说明。
[odefile.m]
function varargout = odefile(t,y,flag,p1,p2)
% odefile.m 是 MATLAB 提供的模板文件。用户可根据需要，使用任何其他文件名。
% varargout 是 MATLAB 专门设计的 " 变长度 " 输出宗量。它由元胞数组构成，因此
% 可适应任意多个、任意形式的输出宗量。请参看第 8.5.2 节。
% t 自变量。不管在下面指令中是否出现，自变量 t 必须作为第一输入宗量。
% y 一阶微分方程组的列向量形式函数。它必须作为第二输入宗量。
% flag 切换变量。它必须处在第三输入宗量位置。用户无须也不要对它直接赋
% 值；它的赋值由微分方程解算指令（ solver ），自动产生。
% p1, p2 是被传递的参数。这里，作为示意，仅列出两个。
% 根据需要，传递参数的数目不受限制。

% 以下是 switch-case 构成的多向选择控制。应注意：
% （ A ）各情况的情况表达式是 MATLAB 规定的，不得任意改变。
% （ B ）各情况所执行的任务也是规定的，不得任意改变。
% （ C ）各情况内，变长输出宗量元胞数是规定的，不得任意改变。
% （ D ）各情况内，（等式右边的）函数名称可以改变，但相应子函数名称要一致。
% （ E ）各情况内，（等式右边的）函数中所包含的 t , y 是必须的，
% 不管它们是否以显式出现相应子函数体中。
% （ F ）各情况内，（等式右边的）函数中所包含的 p1, p2 是传递参数的示意性表示。
% 具体参数视需要而定，只要该参数已经传入 odefile.m 函数内存。

switch flag
case '' % 规定空字符串 '' 情况：专管一阶导数 dy/dt = f(t,y) 的计算
varargout{1} = f(t,y,p1,p2);
case 'init' % 规定 'init' 情况：专管三个宗量 [tspan , y0 , options] 的设置。
[varargout{1:3}] = init(p1,p2);
case 'jacobian' % 规定 'jacobian' 情况：专管计算解析的 Jacobian 矩阵 df/dy 。
varargout{1} = jacobian(t,y,p1,p2);
case 'jpattern' % 规定 'jpattern' 情况：专管计算稀疏的数值 Jacobian 矩阵 df/dy 。
varargout{1} = jpattern(t,y,p1,p2);
case 'mass ' % 规定 'mass' 情况：专管计算质量矩阵（ mass matrix ）。
varargout{1} = mass(t,y,p1,p2);
case 'events' % 规定 'events' 情况：专管事件定义和判断
[varargout{1:3}] = events(t,y,p1,p2);
otherwise
error([ 'Unknown flag ''' flag '''.' ]);
end

% 以下是 odefile.m 的子函数，它们分别与 " 父 " 函数中各情况对应。
%---------------------------------------------------------------------------
function dydt = f(t,y,p1,p2)
% 空字符串 '' 情况内调用的子函数：专管一阶导数 dy/dt 的计算
% dydt 一阶导数。该变量名可由用户按需要取名。
% 下面函数内容，由用户自己编写，最后产生输出 dydt 。
dydt = 〈由用户编写〉
% ---------------------------------------------------------------------------
function [tspan,y0,options] = init(p1,p2)
% 'init' 情况内调用的子函数：专管三个宗量 [tspan , y0 , options] 的设置。
% 三个输出宗量的名称可由用户自起，但各位置上宗量的性质不能改变。
% 下面函数内容，由用户自己编写，最后产生输出三个输出。
tspan = 〈积分的时间区段或时间点集〉
y0 = 〈初值〉
options = 〈或用 odeset 设置算法参数，或用空阵符 [ ] 表示缺省设置〉
% ---------------------------------------------------------------------------

function dfdy = jacobian(t,y,p1,p2)
% 'jacobian' 情况调用的子函数：专管计算解析的 Jacobian 矩阵 df/dy 。
% 输出宗量的名称可由用户自起。
% 下面函数内容，由用户自己编写，最后产生输出输出 dfdy 。
dfdy = 〈 Jacobian 矩阵〉
% ---------------------------------------------------------------------------
function S = jpattern(t,y,p1,p2)
% 'jpattern' 情况调用的子函数：专管计算解析的 Jacobian 矩阵 df/dy 。
% 输出宗量的名称可由用户自起。
% 下面函数内容，由用户自己编写，最后产生输出 S 。
S = 〈稀疏形式的数值 Jacobian 矩阵〉
% ---------------------------------------------------------------------------

function M = mass(t,y,p1,p2)
% 'mass' 情况调用的子函数：专管计算质量矩阵（ mass matrix ）。解刚性方程时用。
% 输出宗量的名称可由用户自起。
% 下面函数内容，由用户自己编写，最后产生输出 M 。
M = 〈质量矩阵〉
% --------------------------------------------------------------------------
function [value,isterminal,direction] = events(t,y,p1,p2)
% 'events' 情况内调用的子函数：专管三个宗量 [value,isterminal,direction] 的设置。
% 三个输出宗量的名称可由用户自起，但各位置上宗量的性质不能改变。
% 下面函数内容，由用户自己编写，最后产生输出三个输出。
value = 〈各元素是标量事件函数，即函数穿越 0 点为事件〉
direction = 〈各元素定义事件函数穿越方式： 1 为向正穿越 0 ； -1 为向负穿越 0 ； 0 不管方向〉
isterminal = 〈各元素定义事件发生后计算是否继续： 1 为终止计算； 0 为继续计算〉
% --------------------------------------------------------------------------

（2） ODE 模板的使用方法

5.14.4 关于解算指令选项options的属性设置
5.14.4.2 options属性处理和输出函数使用演示

【 * 例 5.14.4 .2-1 】仍以卫星轨道问题为例（原题见例 5.14.2.1-1, 5.14.2.1-2 , 5.14.2.2-1 ）。本例演示如何通过对 options 域的直接设置，借助微分方程解算输出指令，表现解算的中间结果。具体目标是：画出解向量中由构成的相平面。相平面的绘制是在微分方程解算中间逐步完成的。

（1）编写 ODE 文件 DYDt4.m （在 DYDt3.m 基础上删改而成）
[DYDt4.m]
function varargout=DYDt4(t,Y,flag,G,ME,tspan,Y0)
switch flag
case ''
varargout{1} = f(t,Y,G,ME);
case 'init'
[varargout{1:3}] = fi(tspan,Y0);
otherwise
error([ 'Unknown flag ''' flag '''.' ]);
end
% ------------------------------------------------------------------
function Yd = f(t,Y,G,ME)
X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V; -G*ME*X/r^3];
% ------------------------------------------------------------------
function [ts,y0,options] = fi(tspan,Y0)
ts=tspan;y0 = Y0;
% 采用向域直接赋值法，设置 options 属性。以供与 odeset 使用方法对照。
options.RelTol=1e-5;options.AbsTol=1e-4;
options.OutputFcn= 'odephas2' ; % 在积分进程中，绘制相平面图。
options.OutputSel=[1 3]; % 解向量的第 1 、 3 分量分别为相平面图的横、纵坐标量。

（2）编写脚本文件 odeexp4.m
[odeexp4.m]
%odeexp4.m
G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;
tspan=[t0,tf];Y0=[x0;0;0;vy0];
% 以下四条指令，为相平面图预置一个坐标轴范围
clf,set(gca, 'xlim' ,[-5 25]*1e7, 'ylim' ,[-3 3]*1e3); % 设置适当坐标范围 <5>
box on % 框形坐标
hold on; % 使中间结果绘在同一幅图上
ode45( 'DYDt4' ,[],[],[],G,ME,tspan,Y0);hold off

（3）在 MATLAB 指令窗中运行以下指令，即可在图形窗中见到相平面图的逐步绘制。
shg,odeexp4