ode15i

Solve fully implicit differential equations — variable order method

collapse all in page

Syntax

[t,y] = ode15i(odefun,tspan,y0,yp0)

[t,y] = ode15i(odefun,tspan,y0,yp0,options)

[t,y,te,ye,ie] = ode15i(odefun,tspan,y0,yp0,options)

sol = ode15i(___)

Description

example

[t,y] = ode15i(odefun,tspan,y0,yp0), wheretspan = [t0 tf], integrates the system of differential equations $f (t, y, y') = 0$ fromt0totfwith initial conditionsy0andyp0. Each row in the solution arrayycorresponds to a value returned in column vectort.

example

[t,y] = ode15i(odefun,tspan,y0,yp0,options)also uses the integration settings defined byoptions, which is an argument created using theodesetfunction. For example, use theAbsTolandRelToloptions to specify absolute and relative error tolerances, or theJacobianoption to provide the Jacobian matrix.

[t,y,te,ye,ie] = ode15i(odefun,tspan,y0,yp0,options)additionally finds where functions of(t,y,y'), called event functions, are zero. In the output,teis the time of the event,yeis the solution at the time of the event, andieis the index of the triggered event.

For each event function, specify whether the integration is to terminate at a zero and whether the direction of the zero crossing matters. Do this by setting the'Events'property to a function, such asmyEventFcnor@myEventFcn, and creating a corresponding function: [value,isterminal,direction] =myEventFcn(t,y,yp). For more information, seeODE Event Location.

sol= ode15i(___)returns a structure that you can use withdevalto evaluate the solution at any point on the interval[t0 tf]. You can use any of the input argument combinations in previous syntaxes.

Examples

collapse all

解决Weissinger隐含的颂歌

Open Live Script

Calculate consistent initial conditions and solve an implicit ODE withode15i.

Weissinger's equation is

${ty}^{2} {(y^{'})}^{3} - y^{3} {(y^{'})}^{2} + t (t^{2} + 1) y^{'} - t^{2} y = 0$ .

Since the equation is in the generic form $f (t, y, y^{'}) = 0$ , you can use theode15ifunction to solve the implicit differential equation.

Code Equation

To code the equation in a form suitable forode15i, you need to write a function with inputs for $t$ , $y$ , and $y^{'}$ that returns the residual value of the equation. The function@weissingerencodes this equation. View the function file.

typeweissinger

function res = weissinger(t,y,yp) %WEISSINGER Evaluate the residual of the Weissinger implicit ODE % % See also ODE15I. % Jacek Kierzenka and Lawrence F. Shampine % Copyright 1984-2014 The MathWorks, Inc. res = t*y^2 * yp^3 - y^3 * yp^2 + t*(t^2 + 1)*yp - t^2 * y;

Calculate Consistent Initial Conditions

Theode15isolver requiresconsistentinitial conditions, that is, the initial conditions supplied to the solver must satisfy

$f (t_{0}, y, y^{'}) = 0$ .

Since it is possible to supply inconsistent initial conditions, andode15idoes not check for consistency, it is recommended that you use the helper functiondecicto compute such conditions.decicholds some specified variables fixed and computes consistent initial values for the unfixed variables.

In this case, fix the initial value $y (t_{0}) = \sqrt{\frac{3}{2}}$ and letdeciccompute a consistent initial value for the derivative $y^{'} (t_{0})$ , starting from an initial guess of $y^{'} (t_{0}) = 0$ .

t0 = 1; y0 = sqrt(3/2); yp0 = 0; [y0,yp0] = decic(@weissinger,t0,y0,1,yp0,0)

y0 = 1.2247

yp0 = 0.8165

Solve Equation

Use the consistent initial conditions returned bydecicwithode15ito solve the ODE over the time interval $[1 10]$ .

[t,y] = ode15i(@weissinger,[1 10],y0,yp0);

Plot Results

The exact solution of this ODE is

$y (t) = \sqrt{t^{2} + \frac{1}{2}}$ .

Plot the numerical solutionycomputed byode15iagainst the analytical solutionytrue.

ytrue = sqrt(t.^2 + 0.5); plot(t,y,'*',t,ytrue,'-o') legend('ode15i','exact')

Figure contains an axes object. The axes object contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent ode15i, exact.

Solve Robertson Problem as Implicit Differential Algebraic Equations (DAEs)

Open Script

This example reformulates a system of ODEs as a fully implicit system of differential algebraic equations (DAEs). The Robertson problem coded byhb1ode.mis a classic test problem for programs that solve stiff ODEs. The system of equations is

$\begin{array}{cl} y'_1 &= -0.04y_1 + 10^4 y_2y_3\\ y'_2 &= 0.04y_1 - 10^4 4y_2y_3-(3 \times 10^7)y_2^2\\ y'_3 &= (3 \times 10^7)y_2^2.\end{array}$

hb1odesolves this system of ODEs to steady state with the initial conditions $y_1 = 1$ , $y_2 = 0$ , and $y_3 = 0$ . But the equations also satisfy a linear conservation law,

$$y'_1 + y'_2 + y'_3 = 0.$$

In terms of the solution and initial conditions, the conservation law is

$$y_1 + y_2 + y_3 = 1.$$

The problem can be rewritten as a system of DAEs by using the conservation law to determine the state of $y_3$ . This reformulates the problem as the implicit DAE system

$\begin{array}{cl} 0 &= y'_1 +0.04y_1 - 10^4 y_2y_3\\ 0 &= y'_2 -0.04y_1 + 10^4 y_2y_3+(3 \times 10^7)y_2^2\\ 0 &= y_1 + y_2 + y_3 - 1.\end{array}$

The functionrobertsidaeencodes this DAE system.

functionres = robertsidae(t,y,yp) res = [yp(1) + 0.04*y(1) - 1e4*y(2)*y(3); yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2; y(1) + y(2) + y(3) - 1];

The full example code for this formulation of the Robertson problem is available inihb1dae.m.

Set the error tolerances and the value of $\partial f / \partial y'$ .

options = odeset('RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6],...'Jacobian',{[],[1 0 0; 0 1 0; 0 0 0]});

Usedecicto compute consistent initial conditions from guesses. Fix the first two components ofy0to get the same consistent initial conditions as found byode15sinhb1dae.m, which formulates this problem as a semi-explicit DAE system.

y0 = [1; 0; 1e-3]; yp0 = [0; 0; 0]; [y0,yp0] = decic(@robertsidae,0,y0,[1 1 0],yp0,[],options);

Solve the system of DAEs usingode15i.

tspan = [0 4 * logspace (6,6)];[t、y] = ode15i (@robertsidae,tspan,y0,yp0,options);

Plot the solution components. Since the second solution component is small relative to the others, multiply it by1e4before plotting.

y(:,2) = 1e4*y(:,2); semilogx(t,y) ylabel('1e4 * y(:,2)') title('Robertson DAE problem with a Conservation Law, solved by ODE15I')

Input Arguments

collapse all

`odefun`—Functions to solve
function handle

Functions to solve, specified as a function handle that defines the functions to be integrated.

The functionf = odefun(t,y,yp), for a scalartand column vectorsyandyp, must return a column vectorfof data typesingleordoublethat corresponds to $f (t, y, y')$ .odefunmust accept the three inputs fort,y, andypeven if one of the inputs is not used in the function.

For example, to solve $y' - y = 0$ , use this function.

functionf = odefun(t,y,yp) f = yp - y;end

For a system of equations, the output ofodefunis a vector. Each equation becomes an element in the solution vector. For example, to solve

$\begin{array}{l} y'_{1} - y_{2} = 0 \\ y'_{2} + 1 = 0, \end{array}$

use this function.

functiondy = odefun(t,y,yp) dy = zeros(2,1); dy(1) = yp(1)-y(2); dy(2) = yp(2)+1;end

For information on how to provide additional parameters to the functionodefun, seeParameterizing Functions.

Example:@myFcn

Data Types:function_handle

`tspan`—Interval of integration
vector

Interval of integration, specified as a vector. At minimum,tspanmust be a two element vector[t0 tf]specifying the initial and final times. To obtain solutions at specific times betweent0andtf, use a longer vector of the form[t0,t1,t2,...,tf]. The elements intspanmust be all increasing or all decreasing.

The solver imposes the initial conditions given byy0at the initial timetspan(1), then integrates fromtspan(1)totspan(结束):

Iftspanhas two elements,[t0 tf],那么解算器返回解决方案评估t each internal integration step within the interval.
Iftspanhas more than two elements[t0,t1,t2,...,tf],那么解算器返回解决方案评估t the given points. However, the solver does not step precisely to each point specified intspan. Instead, the solver uses its own internal steps to compute the solution, then evaluates the solution at the requested points intspan. The solutions produced at the specified points are of the same order of accuracy as the solutions computed at each internal step.
Specifying several intermediate points has little effect on the efficiency of computation, but for large systems it can affect memory management.

The values oftspanare used by the solver to calculate suitable values forInitialStepandMaxStep:

Iftspancontains several intermediate points[t0,t1,t2,...,tf], then the specified points give an indication of the scale for the problem, which can affect the value ofInitialStepused by the solver. Therefore, the solution obtained by the solver might be different depending on whether you specifytspanas a two-element vector or as a vector with intermediate points.
The initial and final values intspanare used to calculate the maximum step sizeMaxStep. Therefore, changing the initial or final values intspancould lead to the solver using a different step sequence, which might change the solution.

Example:[1 10]

Example:[1 3 5 7 9 10]

Data Types:single|double

`y0`—Initial conditions fory
vector

Initial conditions fory, specified as a vector.y0must be the same length as the vector output ofodefun, so thaty0contains an initial condition for each equation defined inodefun.

The initial conditions fory0andyp0must be consistent, meaning that $f (t_{0}, y_{0}, y'_{0}) = 0$ . Use thedecicfunction to compute consistent initial conditions close to guessed values.

Data Types:single|double

`yp0`—Initial conditions fory’
column vector

Initial conditions fory’, specified as a column vector.yp0must be the same length as the vector output ofodefun, so thatyp0contains an initial condition for each variable defined inodefun.

The initial conditions fory0andyp0must be consistent, meaning that $f (t_{0}, y_{0}, y'_{0}) = 0$ . Use thedecicfunction to compute consistent initial conditions close to guessed values.

Data Types:single|double

`options`—Option structure
structure array

Option structure, specified as a structure array. Use theodesetfunction to create or modify the option structure.

SeeSummary of ODE Optionsfor a list of which options are compatible with each ODE solver.

Example:options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot)specifies a relative error tolerance of1e-5, turns on the display of solver statistics, and specifies the output function@odeplotto plot the solution as it is computed.

Data Types:struct

Output Arguments

collapse all

`t`— Evaluation points
column vector

Evaluation points, returned as a column vector.

Iftspancontains two elements[t0 tf], thentcontains the internal evaluation points used to perform the integration.
Iftspancontains more than two elements, thentis the same astspan.

`y`— Solutions
array

Solutions, returned as an array. Each row inycorresponds to the solution at the value returned in the corresponding row oft.

`te`— Time of events
column vector

Time of events, returned as a column vector. The event times intecorrespond to the solutions returned inye, andiespecifies which event occurred.

`ye`— Solution at time of events
array

Solution at time of events, returned as an array. The event times intecorrespond to the solutions returned inye, andiespecifies which event occurred.

`ie`— Index of triggered event function
column vector

Index of triggered event function, returned as a column vector. The event times intecorrespond to the solutions returned inye, andiespecifies which event occurred.

`sol`— Structure for evaluation
structure array

Structure for evaluation, returned as a structure array. Use this structure with thedevalfunction to evaluate the solution at any point in the interval[t0 tf]. Thesolstructure array always includes these fields:

Structure Field	Description
`sol.x`	Row vector of the steps chosen by the solver.
`sol.y`	Solutions. Each column`sol.y(:,i)`contains the solution at time`sol.x(i)`.
`sol.solver`	Solver name.

Additionally, if you specify theEventsoption ofodesetand events are detected, thensolalso includes these fields:

Structure Field	Description
`sol.xe`	Points when events occurred.`sol.xe(end)`contains the exact point of a terminal event, if any.
`sol.ye`	Solutions that correspond to events in`sol.xe`.
`sol.ie`	Indices into the vector returned by the function specified in the`Events`option. The values indicate which event the solver detected.

Tips

Providing the Jacobian matrix toode15iis critical for reliability and efficiency. Alternatively, if the system is large and sparse, then providing the Jacobian sparsity pattern also assists the solver. In either case, useodesetto pass in the matrices using theJacobianorJPatternoptions.

Algorithms

ode15iis a variable-step, variable-order (VSVO) solver based on the backward differentiation formulas (BDFs) of orders 1 to 5.ode15iis designed to be used with fully implicit differential equations and index-1 differential algebraic equations (DAEs). The helper functiondeciccomputes consistent initial conditions that are suitable to be used withode15i[1].

References

[1] Lawrence F. Shampine, “Solving 0 = F(t, y(t), y′(t)) in MATLAB,”Journal of Numerical Mathematics, Vol.10, No.4, 2002, pp. 291-310.

Version History

Introduced before R2006a

ode15i

Syntax

Description

Examples

解决Weissinger隐含的颂歌

Solve Robertson Problem as Implicit Differential Algebraic Equations (DAEs)

Input Arguments

`odefun`—Functions to solve
function handle

`tspan`—Interval of integration
vector

`y0`—Initial conditions fory
vector

`yp0`—Initial conditions fory’
column vector

`options`—Option structure
structure array

Output Arguments

`t`— Evaluation points
column vector

`y`— Solutions
array

`te`— Time of events
column vector

`ye`— Solution at time of events
array

`ie`— Index of triggered event function
column vector

`sol`— Structure for evaluation
structure array

Tips

Algorithms

References

Version History

See Also

Topics

ode15i

Syntax

Description

Examples

解决Weissinger隐含的颂歌

Solve Robertson Problem as Implicit Differential Algebraic Equations (DAEs)

Input Arguments

odefun—Functions to solvefunction handle

tspan—Interval of integrationvector

y0—Initial conditions foryvector

yp0—Initial conditions fory’column vector

options—Option structurestructure array

Output Arguments

t— Evaluation pointscolumn vector

y— Solutionsarray

te— Time of eventscolumn vector

ye— Solution at time of eventsarray

ie— Index of triggered event functioncolumn vector

sol— Structure for evaluationstructure array

Tips

Algorithms

References

Version History

See Also

Topics

`odefun`—Functions to solve
function handle

`tspan`—Interval of integration
vector

`y0`—Initial conditions fory
vector

`yp0`—Initial conditions fory’
column vector

`options`—Option structure
structure array

`t`— Evaluation points
column vector

`y`— Solutions
array

`te`— Time of events
column vector

`ye`— Solution at time of events
array

`ie`— Index of triggered event function
column vector

`sol`— Structure for evaluation
structure array