decic

Compute consistent initial conditions for ode15i

Syntax

[y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0)

[y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0,options)

[y0_new,yp0_new,resnrm] = decic(___)

Description

[y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0) uses y0 and yp0 as guesses for the initial conditions of the fully implicit function odefun, holds the components specified by fixed_y0 and fixed_yp0 as fixed, then computes values for the nonfixed components. The result is a complete set of consistent initial conditions. The new values yo_new and yp0_new satisfy odefun(t0,y0_new,yp0_new) = 0 and are suitable to be used as initial conditions with ode15i.

[y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0,options) also uses the options structure options to specify values for AbsTol and RelTol. Create the options structure using odeset.

[y0_new,yp0_new,resnrm] = decic(___) returns the norm of odefun(t0,y0_new,yp0_new) as resnrm. If the norm seems unduly large, then use options to decrease the relative error tolerance RelTol, which has a default value of 1e-3.

Examples

Compute Consistent Initial Conditions for Implicit Equations

Consider the implicit system of equations

$$\begin{array}{cc}0& =2y_1^{\prime }-y_2\\ 0& =y_1+y_2\end{array}$$

These equations are straightforward enough that it is simple to read off consistent initial conditions for the variables. For example, if you fix $$y_1=1$$, then $$y_2=-1$$ according to the second equation and $$y_1^{\prime }=-1/2$$ according to the first equation. Since these values of $$y_1$$, $$y_1^{\prime }$$, and $$y_2$$ satisfy the equations, they are consistent.

Confirm these values by using decic to compute consistent initial conditions for the equations, fixing the value $$y_1=1$$. Use guesses of y0 = [1 0] and yp0 = [0 0], which do not satisfy the equations and are thus inconsistent.

odefun = @(t,y,yp) [2*yp(1)-y(2); y(1)+y(2)];
t0 = 0;
y0 = [1 0];
yp0 = [0 0];
[y0,yp0] = decic(odefun,t0,y0,[1 0],yp0,[])

y0 = 2×1

     1
    -1

yp0 = 2×1

   -0.5000
         0

Solve Weissinger Implicit ODE

Calculate consistent initial conditions and solve an implicit ODE with ode15i.

Weissinger's equation is

${\mathit{t}\mathit{y}}^2{\left({\mathit{y}}^{\prime }\right)}^3-{\mathit{y}}^3{\left({\mathit{y}}^{\prime }\right)}^2+\mathit{t}\left({\mathit{t}}^2+1\right){\mathit{y}}^{\prime }-{\mathit{t}}^2\mathit{y}=0$.

Since the equation is in the generic form $\mathit{f}\left(\mathit{t},\mathit{y},{\mathit{y}}^{\prime }\right)=0$, you can use the ode15i function to solve the implicit differential equation.

Code Equation

To code the equation in a form suitable for ode15i, you need to write a function with inputs for $\mathit{t}$, $\mathit{y}$, and ${\mathit{y}}^{\prime }$ that returns the residual value of the equation. The function @weissinger encodes this equation. View the function file.

type weissinger

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

The ode15i solver requires consistent initial conditions, that is, the initial conditions supplied to the solver must satisfy

$\mathit{f}\left({\mathit{t}}_0,\mathit{y},{\mathit{y}}^{\prime }\right)=0$.

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

In this case, fix the initial value $\mathit{y}\left({\mathit{t}}_0\right)=\sqrt{\frac{3}{2}}$ and let decic compute a consistent initial value for the derivative ${\mathit{y}}^{\prime }\left({\mathit{t}}_0\right)$, starting from an initial guess of ${\mathit{y}}^{\prime }\left({\mathit{t}}_0\right)=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 by decic with ode15i to solve the ODE over the time interval $\left[110\right]$.

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

Plot Results

The exact solution of this ODE is

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

Plot the numerical solution y computed by ode15i against the analytical solution ytrue.

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

Input Arguments

odefun — Functions to solve

function handle

Functions to solve, specified as a function handle that defines the functions to be integrated. odefun represents the system of implicit differential equations that you want to solve using ode15i.

The function f = odefun(t,y,yp), for a scalar t and column vectors y and yp, must return a column vector f of data type single or double that corresponds to $$f\left(t,y,y\text{'}\right)$$. odefun must accept all three input arguments, t, y, and yp even if one of the arguments is not used in the function.

For example, to solve $$y\text{'}-y=0$$, use this function.

function f = odefun(t,y,yp)
f = yp - y;

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

$$\begin{array}{l}y{\text{'}}_1-y_2=0\\ y{\text{'}}_2+1=0,\end{array}$$

use this function.

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

For information on how to provide additional parameters to the function odefun, see Parameterizing Functions.

Example: @myFcn

Data Types: function_handle

t0 — Initial time

scalar

Initial time, specified as a scalar. decic uses the initial time to compute consistent initial conditions that satisfy odefun(t0,y0_new,yp0_new) = 0.

Data Types: single | double

y0 — Initial guesses for y-components

vector

Initial guesses for y-components, specified as a vector. Each element in y0 specifies an initial condition for one dependent variable $$y_n$$ in the system of equations defined by odefun.

Data Types: single | double

fixed_y0 — y-components to hold fixed

vector of 1s and 0s | []

y-components to hold fixed, specified as a vector of 1s and 0s, or as [].

• Set fixed_y0(i) = 1 if no change is permitted in the guess for y0(i).

• Set fixed_y0 = [] if any entry can be changed.

You cannot fix more than length(yp0) components. Depending on the specific problem, it is not always possible to fix certain components of y0 or yp0. It is a best practice not to fix more components than is necessary.

yp0 — Initial guesses for y'-components

vector

Initial guesses for y'-components, specified as a vector. Each element in yp0 specifies an initial condition for one differentiated dependent variable $$y{\text{'}}_n$$ in the system of equations defined by odefun.

Data Types: single | double

fixed_yp0 — y'-components to hold fixed

vector of 1s and 0s | []

y'-components to hold fixed, specified as a vector of 1s and 0s, or as [].

• Set fixed_yp0(i) = 1 if no change is permitted in the guess for yp0(i).

• Set fixed_yp0 = [] if any entry can be changed.

You cannot fix more than length(yp0) components. Depending on the specific problem, it is not always possible to fix certain components of y0 or yp0. It is a best practice not to fix more components than is necessary.

options — Options structure

structure array

Options structure, specified as a structure array. Use the odeset function to create or modify the options structure. The relevant options for use with the decic function are RelTol and AbsTol, which control the error thresholds used to compute the initial conditions.

Example: options = odeset('RelTol',1e-5)

Data Types: struct

Output Arguments

y0_new — Consistent initial conditions for y0

vector

Consistent initial conditions for y0, returned as a vector. If the value of resnrm is small, then yo_new and yp0_new satisfy odefun(t0,y0_new,yp0_new) = 0 and are suitable to be used as initial conditions with ode15i.

yp0_new — Consistent initial conditions for yp0

vector

Consistent initial conditions for yp0, returned as a vector. If the value of resnrm is small, then yo_new and yp0_new satisfy odefun(t0,y0_new,yp0_new) = 0 and are suitable to be used as initial conditions with ode15i.

resnrm — Norm of residual

vector

Norm of residual, returned as a vector. resnrm is the norm of odefun(t0,y0_new,yp0_new).

• A small value of resnrm indicates that decic successfully computed consistent initial conditions that satisfy odefun(t0,y0_new,yp0_new) = 0.

• If the value of resnrm is large, try adjusting the error thresholds RelTol and AbsTol using the options input.

Tips

• The ihb1dae and iburgersode example files use decic to compute consistent initial conditions before solving with ode15i. Type edit ihb1dae or edit iburgersode to view the code.

• You can additionally use decic to compute consistent initial conditions for DAEs solved by ode15s or ode23t. To do this, follow these steps.

  1. Rewrite the system of equations in fully implicit form f(t,y,y') = 0.

  2. Call decic to compute consistent initial conditions for the equations.

  3. Specify y0_new as the initial condition in the call to the solver, and specify yp_new as the value of the InitialSlope option of odeset.

Version History

Introduced before R2006a