Solve DAEs Using Mass Matrix Solvers
Solve differential algebraic equations by using one of the mass matrix solvers available in MATLAB®. To use this workflow, first complete steps 1, 2, and 3 from Solve Differential Algebraic Equations (DAEs). Then, use a mass matrix solver instead of ode15i
.
This example demonstrates the use of ode15s
or ode23t
. For details on the other solvers, see Choose an ODE Solver and adapt the workflow on this page.
Step 1. Convert DAEs to Function Handles
The code from steps 1 to 3 of Solve Differential Algebraic Equations (DAEs) is:
syms x(t) y(t) T(t) m r g eqn1 = m*diff(x(t), 2) == T(t)/r*x(t); eqn2 = m*diff(y(t), 2) == T(t)/r*y(t) - m*g; eqn3 = x(t)^2 + y(t)^2 == r^2; eqns = [eqn1 eqn2 eqn3]; vars = [x(t); y(t); T(t)]; origVars = length(vars); [eqns,vars] = reduceDifferentialOrder(eqns,vars); [DAEs,DAEvars] = reduceDAEIndex(eqns,vars); [DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars);
From the output of reduceDAEIndex
and reduceRedundancies
, you have a vector of equations DAEs
and a vector of variables DAEvars
. To use ode15s
or ode23t
, you need two function handles: one representing the mass matrix of a DAE system, and the other representing the right sides of the mass matrix equations. These function handles form the equivalent mass matrix representation of the ODE system where .
Find these function handles by using massMatrixForm
to get the mass matrix M
and the right sides F
.
[M,f] = massMatrixForm(DAEs,DAEvars)
M =
f =
The equations in DAEs
can contain symbolic parameters that are not specified in the vector of variables DAEvars
. Find these parameters by using setdiff
on the output of symvar
from DAEs
and DAEvars
.
pDAEs = symvar(DAEs); pDAEvars = symvar(DAEvars); extraParams = setdiff(pDAEs,pDAEvars)
extraParams =
The mass matrix M
does not have these extra parameters. Therefore, convert M
directly to a function handle by using odeFunction
.
M = odeFunction(M,DAEvars);
Convert f
to a function handle. Specify the extra parameters as additional inputs to odeFunction
.
f = odeFunction(f,DAEvars,g,m,r);
The rest of the workflow is purely numerical. Set parameter values and create the function handle.
gVal = 9.81; mVal = 1; rVal = 1; F = @(t,Y) f(t,Y,gVal,mVal,rVal);
Step 2. Find Initial Conditions
The solvers require initial values for all variables in the function handle. Find initial values that satisfy the equations by using the MATLAB decic
function. The decic
accepts guesses (which might not satisfy the equations) for the initial conditions, and tries to find satisfactory initial conditions using those guesses. decic
can fail, in which case you must manually supply consistent initial values for your problem.
First, check the variables in DAEvars
.
DAEvars
DAEvars =
Here, Dxt(t)
is the first derivative of x(t)
, Dytt(t)
is the second derivative of y(t)
, and so on. There are 7 variables in a 7
-by-1
vector. Thus, guesses for initial values of variables and their derivatives must also be 7
-by-1
vectors.
Assume the initial angular displacement of the pendulum is 30° or pi/6
, and the origin of the coordinates is at the suspension point of the pendulum. Given that we used a radius rVal
of 1
, the initial horizontal position x(t)
is rVal*sin(pi/6)
. The initial vertical position y(t)
is -rVal*cos(pi/6)
. Specify these initial values of the variables in the vector y0est
.
Arbitrarily set the initial values of the remaining variables and their derivatives to 0
. These are not good guesses. However, they suffice for our problem. In your problem, if decic
errors, then provide better guesses and refer to the decic
page.
y0est = [rVal*sin(pi/6); -rVal*cos(pi/6); 0; 0; 0; 0; 0]; yp0est = zeros(7,1);
Create an option set that contains the mass matrix M
and initial guesses yp0est
, and specifies numerical tolerances for the numerical search.
opt = odeset(Mass=M,InitialSlope=yp0est, ...
RelTol=1e-7,AbsTol=1e-7);
Find consistent initial values for the variables and their derivatives by using the MATLAB decic
function. The first argument of decic
must be a function handle describing the DAE as f(t,y,yp) = f(t,y,y') = 0
. In terms of M
and F
, this means f(t,y,yp) = M(t,y)*yp - F(t,y)
.
implicitDAE = @(t,y,yp) M(t,y)*yp - F(t,y); [y0,yp0] = decic(implicitDAE,0,y0est,[],yp0est,[],opt)
y0 = 7×1
0.4771
-0.8788
-8.6214
0
-0.0000
-2.2333
-4.1135
yp0 = 7×1
0
-0.0000
0
0
-2.2333
0
0
Now create an option set that contains the mass matrix M
of the system and the vector yp0
of consistent initial values for the derivatives. You will use this option set when solving the system.
opt = odeset(opt,InitialSlope=yp0);
Step 3. Solve DAE System
Solve the system integrating over the time span 0
≤ t
≤ 0.5
. Add the grid lines and the legend to the plot. The code uses ode15s
. Instead, you can also use ode23t
by replacing ode15s
with ode23t
.
[tSol,ySol] = ode15s(F,[0, 0.5],y0,opt); plot(tSol,ySol(:,1:origVars),"-o") for k = 1:origVars S{k} = char(DAEvars(k)); end legend(S,Location="best") grid on
Solve the system for different parameter values by setting the new value and regenerating the function handle and initial conditions.
Set rVal
to 2
and regenerate the function handle and initial conditions.
rVal = 2; F = @(t,Y) f(t,Y,gVal,mVal,rVal); y0est = [rVal*sin(pi/6); -rVal*cos(pi/6); 0; 0; 0; 0; 0]; implicitDAE = @(t,y,yp) M(t,y)*yp - F(t,y); [y0,yp0] = decic(implicitDAE,0,y0est,[],yp0est,[],opt); opt = odeset(opt,InitialSlope=yp0);
Solve the system for the new parameter value.
[tSol,ySol] = ode15s(F,[0, 0.5],y0,opt); plot(tSol,ySol(:,1:origVars),"-o") for k = 1:origVars S{k} = char(DAEvars(k)); end legend(S,Location="best") grid on