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.

`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)

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);

`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)

y0 = 0.5000 -0.8660 -8.4957 0 0 yp0 = 0 0 0 -4.2479 -2.4525

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);

`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*cos(pi/6); -r*sin(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

`daeFunction`

| `decic`

| `findDecoupledBlocks`

| `incidenceMatrix`

| `isLowIndexDAE`

| `massMatrixForm`

| `odeFunction`

| `reduceDAEIndex`

| `reduceDAEToODE`

| `reduceDifferentialOrder`

| `reduceRedundancies`