Solves a set of first-order, differential algebraic equations using the Backward Differentiation Formula.  At each time step, it solves the nonlinear system by Newton iteration.

Syntax

x = DAE11a(func, init, deriv_init, time, reltol, abstol, userdata)

x = DAE11a(func, jacobian, init, deriv_init, start, end, points, reltol, abstol, userdata)

Arguments

Name

Description

func

The name of the function that describes the systems of equations.  Must be a string.  The function must return the residuals of the implicit DAE equations.  See the Comments section.

Jacobian

(optional)

The name of the function that computes the modified Jacobian matrix of the model.  Must be a string. The matrix must be N by N, where N is the number of state variables.  See Example 2 for an explanation of the modified Jacobian.  Also, see the Comments section.

init

A vector of the initial conditions for the state variables of the system.  Its length must equal the number of equations in the system.  The order of the variables will be followed in the system function. See the Comments section.

deriv_init

A vector of the initial conditions for the first derivatives of the state variables of the system.  Its length must equal the length of init.  See the Comments section.

time

A vector of time stamps where the solutions are to be returned.  Must be positive monotonic.

start

A non-negative scalar specifying the start time of the solution.

end

A non-negative scalar specifying the end time of the solution.  Must be greater than start.

points

A positive scalar specifying the number of points between start and end where the solutions are to be returned.

reltol

(optional)

Relative tolerance of convergence.  This applies to all states.  Must be a positive scalar.  Defaults to 1.0e-3.

abstol

(optional)

Absolute tolerances below which the solutions are unimportant.  This is used when a solution approaches zero.  If it is a vector with the same length as init, then each state has its own tolerance.  If it is a one element vector, then the value applies to each state.  If it is an empty vector or omitted, the default is 1.0e-6 for each state.  Every tolerance value must be positive.

userdata (optional)

Data matrix that is passed to the user defined function at each iteration.  This can be used to supply data that are defined outside the function.  It can be omitted or empty.

Output

Name

Description

x

A matrix of the solutions at each time stamp.  Each column contains the values of the state variables at a requested time.

Example 1

Solve the van der Pol equation for a time span of 0 to 10 seconds with increments of 0.1 second.

Syntax

// First define the van der Pol system in a function

function vdp(t, y, yp, d)

{

 // t is a scalar, the independent variable

 // y is a 2 element state vector

     // yp = dy/dt

     // d is the user data (not used)

 mu = 1;

 // Equation for first residual

 r1 = y(2) - yp(1)

 // Equation for second residual

 r2 = mu * (1.0 - y(1)^2) * y(2) - y(1) - yp(2)

 r = [r1; r2]    // return the residuals vector

 return r

}

yi = [2; 0]

ypi = [0; -2]

t = [0:0.1:10]; // time span

x = DAE11a("vdp", yi, ypi, t, 0.001, [0.001,0.001]);

PlotLine(t, x(:,1));

PlotLine(t, x(:,2));

Result

x will be a 101x2 matrix, with the first column representing the dependant variable and the second column its derivative in this case.

Example 2

Modify the above example to use the Jacobian matrix.

Syntax

// Define the Van der Pol Jacobian

function jacob(t, y, yp, r, c, d)

{

     // yp = dy/dt

     // r is residuals vector

     // c is a step control scalar, not user supplied

     // A modified Jacobian is required as follows

     // J = dr/dy + c * dr/dyp

 mu = 100

 J = Zeros(2,2)

 J(1, 1) = -c

 J(1, 2) = 1.0

 J(2, 1) = -2.0 * mu * y(1) * y(2) - 1.0

 J(2, 2) = mu * (1.0 - y(1)^2) – c

 return J

}

x = DAE11a('vdp', 'jacob', yi, ypi, t, 0.001,

          [0.001,0.001]);

Result

x will be a 101x2 matrix with the first column representing the dependant variable and the second column its derivative in this case.

Comments

DAE11a is not suitable for stiff problems.

The initial condition vectors can be either rows or a columns. Inside the user function, the state variable vectors will be columns.

If an optional argument is supplied, all other optional arguments up until that one must be supplied also.

The user function must accept at least three inputs. The first one is the time-step, the second and third inputs are column vectors of the state variable and their derivative values at that step, respectively.  Optionally, it can accept a fourth argument which is a matrix that can be used to pass in any additional user data.

The Jacobian function must take the time-step, the column vectors of the state variable and their derivative values at that step, respectively.  It must return the Jacobian matrix.

See Also:

RK45

ODE21

ODE22a

HMath-4010: Solving Ordinary Differential and Differential Algebraic Equations