Solve Semilinear DAE System

This workflow is an alternative workflow to solving semilinear DAEs, used only if reduceDAEIndex failed in the standard workflow with the warning below. For the standard workflow, see Solve Differential Algebraic Equations (DAEs).

Warning: The index of the reduced DAEs is larger...
than 1. [daetools::reduceDAEIndex]

To solve your DAE system complete these steps.

Step 1. Reduce Differential Index with `reduceDAEToODE`

Complete steps 1 and 2 in Solve Differential Algebraic Equations (DAEs) before beginning this step. Then, in step 3 when reduceDAEIndex fails, reduce the differential index using reduceDAEToODE. The advantage of reduceDAEToODE is that it reliably reduces semilinear DAEs to ODEs (DAEs of index 0). However, this function is slower and works only on semilinear DAE systems. reduceDAEToODE can fail if the system is not semilinear.

To reduce the differential index of the DAEs described by eqns and vars, use reduceDAEToODE. To reduce the index, reduceDAEToODE adds new variables and equations to the system. reduceDAEToODE also returns constraints, which are conditions that help find initial values to ensure that the resulting ODEs are equal to the initial DAEs.

[ODEs,constraints] = reduceDAEToODE(eqns,vars)

ODEs =
                             Dxt(t) - diff(x(t), t)
                             Dyt(t) - diff(y(t), t)
                  m*diff(Dxt(t), t) - (T(t)*x(t))/r
          m*diff(Dyt(t), t) - (T(t)*y(t) - g*m*r)/r
        -(4*T(t)*y(t) - 2*g*m*r)*diff(y(t), t) -...
           diff(T(t), t)*(2*x(t)^2 + 2*y(t)^2) -...
                     4*T(t)*x(t)*diff(x(t), t) -...
                  4*m*r*Dxt(t)*diff(Dxt(t), t) -...
                       4*m*r*Dyt(t)*diff(Dyt(t), t)
 
constraints =
 2*g*m*r*y(t) - 2*T(t)*y(t)^2 - 2*m*r*Dxt(t)^2 -...
                     2*m*r*Dyt(t)^2 - 2*T(t)*x(t)^2
                              r^2 - y(t)^2 - x(t)^2
                       2*Dxt(t)*x(t) + 2*Dyt(t)*y(t)

Step 2. ODEs to Function Handles for ode15s and ode23t

From the output of reduceDAEToODE, you have a vector of equations in ODEs and a vector of variables in vars. To use ode15s or ode23t, you need two function handles: one representing the mass matrix of the ODE system, and the other representing the vector containing the right sides of the mass matrix equations. These function handles are the equivalent mass matrix representation of the ODE system where M(t,y(t))y’(t) = f(t,y(t)).

Find these function handles by using massMatrixForm to get the mass matrix massM (M in the equation) and right sides f.

[massM,f] = massMatrixForm(ODEs,vars)

massM =
[           -1,                     0,                     0,             0,             0]
[            0,                    -1,                     0,             0,             0]
[            0,                     0,                     0,             m,             0]
[            0,                     0,                     0,             0,             m]
[ -4*T(t)*x(t), 2*g*m*r - 4*T(t)*y(t), - 2*x(t)^2 - 2*y(t)^2, -4*m*r*Dxt(t), -4*m*r*Dyt(t)]
 
f =
               -Dxt(t)
               -Dyt(t)
         (T(t)*x(t))/r
 (T(t)*y(t) - g*m*r)/r
                     0

The equations in ODEs can contain symbolic parameters that are not specified in the vector of variables vars. Find these parameters by using setdiff on the output of symvar from ODEs and vars.

pODEs = symvar(ODEs);
pvars = symvar(vars);
extraParams = setdiff(pODEs, pvars)

extraParams =
[ g, m, r]

The extra parameters that you need to specify are the mass m, radius r, and gravitational constant g.

Convert massM and f to function handles using odeFunction. Specify the extra symbolic parameters as additional inputs to odeFunction.

massM = odeFunction(massM, vars, m, r, g);
f = odeFunction(f, vars, m, r, g);

The rest of the workflow is purely numerical. Set the parameter values and substitute the parameter values in DAEs and constraints.

m = 1;
r = 1;
g = 9.81;
ODEsNumeric = subs(ODEs);
constraintsNumeric = subs(constraints);

Create the function handle suitable for input to ode15s or ode23s.

M = @(t, Y) massM(t, Y, m, r, g);
F = @(t, Y) f(t, Y, m, r, g);

Step 3. Initial Conditions for `ode15s` and `ode23t`

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 vars.

vars

vars =
   x(t)
   y(t)
   T(t)
 Dxt(t)
 Dyt(t)

Here, Dxt(t) is the first derivative of x(t), and so on. There are 5 variables in a 5-by-1 vector. Therefore, guesses for initial values of variables and their derivatives must also be 5-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 r of 1, the initial horizontal position x(t) is r*sin(pi/6). The initial vertical position y(t) is -r*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 this problem. In your problem, if decic errors, then provide better guesses and refer to the decic page.

y0est = [r*sin(pi/6); -r*cos(pi/6); 0; 0; 0];
yp0est = zeros(5,1);

Create an option set that contains the mass matrix M of the system and specifies numerical tolerances for the numerical search.

opt = odeset('Mass', M, 'RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));

Find initial values consistent with the system of ODEs and with the algebraic constraints by using decic. The parameter [1,0,0,0,1] in this function call fixes the first and the last element in y0est, so that decic does not change them during the numerical search. Here, this fixing is necessary to ensure decic finds satisfactory initial conditions.

[y0, yp0] = decic(ODEsNumeric, vars, constraintsNumeric, 0,...
                y0est, [1,0,0,0,1], yp0est, opt)

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 4. Solve an ODE System with `ode15s` or `ode23t`

Solve the system integrating over the time span 0 ≤ t ≤ 0.5. Add the grid lines and the legend to the plot. Use ode23s by replacing ode15s with ode23s.

[tSol,ySol] = ode15s(F, [0, 0.5], y0, opt);
plot(tSol,ySol(:,1:origVars),'-o')

for k = 1:origVars
  S{k} = char(vars(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 r to 2 and repeat the steps.

r = 2;

ODEsNumeric = subs(ODEs);
constraintsNumeric = subs(constraints);
M = @(t, Y) massM(t, Y, m, r, g);
F = @(t, Y) f(t, Y, m, r, g);

y0est = [r*sin(pi/6); -r*cos(pi/6); 0; 0; 0];
opt = odeset('Mass', M, 'RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));
[y0, yp0] = decic(ODEsNumeric, vars, constraintsNumeric, 0,...
		 y0est, [1,0,0,0,1], 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(vars(k));
end

legend(S, 'Location', 'Best')
grid on