ode15i Solve fully
implicit differential equations, variable order method.
[TOUT,YOUT] =
ode15i(ODEFUN,TSPAN,Y0,YP0) with TSPAN = [T0
TFINAL]
integrates the system of
differential equations f(t,y,y') = 0 from
time
T0 to TFINAL, with
initial conditions Y0,YP0. Function ode15i solves
ODEs
and DAEs of index 1. The
initial conditions must be "consistent",
meaning
that f(T0,Y0,YP0)=0. Use
the function DECIC to compute consistent
initial
conditions close to
guessed values. ODEFUN is a function handle. For
a
scalar T and vectors Y
and YP, ODEFUN(T,Y,YP) must return a column vector
corresponding to
f(t,y,y'). Each row in the solution array
YOUT
corresponds to a time
returned in the column vector TOUT. To
obtain
solutions at specific
times T0,T1,...,TFINAL (all increasing or
all
decreasing), use TSPAN =
[T0 T1 ... TFINAL].
[TOUT,YOUT] =
ode15i(ODEFUN,TSPAN,Y0,YP0,OPTIONS) solves as above
with
default integration
properties replaced by values in OPTIONS, an
argument
created with the ODESET
function. See ODESET for details. Commonly
used
options are scalar
relative error tolerance 'RelTol' (1e-3 by
default)
and vector of absolute
error tolerances 'AbsTol' (all components
1e-6
by default).
The Jacobian matrices
df/dy and df/dy' are critical to reliability
and
efficiency.
Use ODESET to set 'Jacobian' to a function handle
FJAC if
[DFDY,DFDYP] =
FJAC(T,Y,YP) returns the matrices df/dy and
df/dy'.
If DFDY = [], df/dy is
approximated by finite differences and
similarly
for DFDYP.
If the 'Jacobian' option is not set (the
default), both
matrices are
approximated by finite differences.
Set 'Vectorized'
{'on','off'} if the ODE function is coded so
that
ODEFUN(T,[Y1 Y2 ...],YP)
returns [ODEFUN(T,Y1,YP) ODEFUN(T,Y2,YP)
...].
Set 'Vectorized'
{'off','on'} if the ODE function is coded so that
ODEFUN(T,Y,[YP1 YP2
...]) returns [ODEFUN(T,Y,YP1) ODEFUN(T,Y,YP2) ...].
If df/dy or df/dy' is a
sparse matrix, set 'JPattern' to the sparsity
patterns, {SPDY,SPDYP}.
A sparsity pattern of df/dy is a sparse
matrix SPDY with
SPDY(i,j) = 1 if component i of f(t,y,yp)
depends on component j
of y, and 0 otherwise. Use SPDY = [] to
indicate that df/dy is a
full matrix. Similarly for df/dy' and
SPDYP. The default value
of 'JPattern' is {[],[]}.
[TOUT,YOUT,TE,YE,IE] =
ode15i(ODEFUN,TSPAN,Y0,YP0,OPTIONS) with the 'Events'
property in OPTIONS set
to a function handle EVENTS, solves as above
while
also finding where
functions of (T,Y,YP), called event functions, are
zero.
For each function you
specify whether the integration is to terminate
at
a zero and whether the
direction of the zero crossing matters. These
are
the three column vectors
returned by EVENTS: [VALUE,ISTERMINAL,DIRECTION] =
EVENTS(T,Y,YP). For the
I-th event function: VALUE(I) is the value of the
function,
ISTERMINAL(I)=1 if the integration is to terminate at a zero
of
this event function and
0 otherwise. DIRECTION(I)=0 if all zeros are to
be computed (the
default), +1 if only zeros where the event function is
increasing, and -1 if
only zeros where the event function is
decreasing. Output TE is
a column vector of times at which events occur.
Rows of YE are the
corresponding solutions, and indices in vector IE
specify which event
occurred.
SOL = ode15i(ODEFUN,[T0
TFINAL],Y0,YP0,...) returns a structure that
can be
used with DEVAL to evaluate the solution or its
first derivative
at any point between T0
and TFINAL. The steps chosen by ode15i are
returned
in a row vector SOL.x.
For each I, the column SOL.y(:,I)
contains
the solution at
SOL.x(I). If events were detected, SOL.xe is a row
vector
of points at which
events occurred. Columns of SOL.ye are the
corresponding
solutions, and indices
in vector SOL.ie specify which event
occurred.
Example
t0 = 1;
y0 = sqrt(3/2);
yp0 = 0;
[y0,yp0] =
decic(@weissinger,t0,y0,1,yp0,0);
uses
a helper function DECIC to hold fixed the initial value for
y(t0)
and
compute a consistent initial value for y'(t0) for the
Weissinger
implicit ODE. The ODE is solved using ode15i and the numerical
solution
is
plotted against the analytical solution
[t,y] = ode15i(@weissinger,[1
10],y0,yp0);
ytrue = sqrt(t.^2 +
0.5);
plot(t,y,t,ytrue,'o');
ode15i()调用格式如下
[y0mod,yp0mod,resnrm] =
decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0, options...)
[T,Y,TE,YE,IE] =
ode15i(odefun,tspan,y0mod,yp0mod,options...)
隐式ODEs不同于一般的显式ODEs，ode15i需要同时给出状态变量及其一阶导数的初值(x0,dx0)，它们不能任意赋值，只能有n个是独立的，其余的需要隐式方程求解(使用decic函数)，否则出现矛盾的初值条件。
输入参数：
odefun：微分方程，注意此时的函数变量有t(自变量)，x(状态变量)，dx(状态变量一阶导数)等三个
t0：自变量初值，必须与tspan(1)相等
y0：状态变量初值，题目给出几个就写几个，没给的自己猜测一个
fixed_y0：与y0的尺寸一样，表示对应位置的初值是给定的还是猜测的，如果给定则1，否则0
yp0：状态变量一阶导数初值，其它同y0
fixed_yp0：同理fixed_y0，不过此时对应的是yp0了
options：其它选项，具体查看帮助
tspan：自变量的范围，其中tspan(1)必须与t0一致
输出参数：
大部分同理ode45.
下面我们以实例说明，对于下面的IDE，求解x
解】
(1)选取状态变量
(2)写出状态变量的一阶微分。只不过此时的x''和y''不再是显式的
(3)书写odefun，注意此时多了一个dx变量，就是x的一阶导数
odefun=@(t,x,dx)[dx(1)-x(2)
dx(2)*sin(x(4))+dx(4)^2+2*x(1)*x(3)-x(1)*dx(2)*x(4)
dx(3)-x(4)
x(1)*dx(2)*dx(4)+cos(dx(4))-3*x(3)*x(2)];
(4)求解dx0
t0=0;%自变量初值
%对于x0和dx0中的题目给出的初值，如实写，没有给出的任意写
x0=[1 0 0 1]';%本题初值x0的都给出了，所以必须是这个
%初值x0中题目给出的，对应位置填1，否则为0
fix_x0=ones(4,1);%本题中x0都给出了，故全为1
dx0=[0 0 1 1]';%本题中初值dx0一个都没有给出，那么全部任意写
%初值dx0中题目给出的，对应位置填1，否则为0
fix_dx0=zeros(4,1);%本题中dx0一个没有给出，故全部为0
[x02,dx02]=decic(odefun,t0,x0,fix_x0,dx0,fix_dx0);
(5)求解微分方程
solution=ode15i(odefun,[0
2],x02,dx02)%x02和dx02是由decic输出参数
虽然我们上面的分析是完全正确，但是好像Matlab就是罢工，郁闷死了
其实对一些简单的ODEs压根没有必要通过decic函数来求解未知的初值，我们可以直接人工算出x02和dx02，根据下面的式子
我们演示下，根据题目x0=[x1,x2,x3,x4]=[1 0 0
1]，由第一个式子可以得出x1’=0，第三个式子有x3’=1，再联立第二四两个式子后有x2’=0，x4’=1，故dx0=[x1’
,x2’,x3’,x4’]=[0 0 1 1]
看了下，隐式ODEs数值求解的确麻烦的些，但是和前面的ode45没有本质的区别，decic函数就是为了求解那些隐式中的状态变量的一阶导数，与我们方法二中的很相似，只不过这里直接调用罢了。
其实不管隐式还是显式ODEs，甚至DAE，都是可以使用ode15i求解，只不过此时将显式当隐式处理，或者说将代数约束看成一个隐式罢了！
参考：http://www.matlabsky.com/thread-529-1-1.html
薛定宇，陈阳泉. <高等应用数学问题的MATLAB求解>
加载中，请稍候......
加载中…
(2015-04-02 19:30:05)