Analyze and Manipulate Differential Algebraic Equations

Define Parameters of the Model

Consider a 2-D physical pendulum, consisting of a mass m attached to the origin by a string of constant length r. Only the gravitational acceleration g = 9.81 m/s^2 acts on the mass. The model consists of second-order differential equation for the position (x(t), y(t)) of the mass with an unknown force F(t) inside the string which serves for keeping the mass on the circle. The force is directed along the string.

syms x(t) y(t) F(t) m g r
eqs = [m*diff(x(t), t, t) == F(t)/r*x(t);
       m*diff(y(t), t, t) == F(t)/r*y(t) - m*g;
       x(t)^2 + y(t)^2 == r^2]

eqs = 
$$\left(\begin{array}{c}m\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }x\left(t\right)=\frac{F\left(t\right)\hspace{0.17em}x\left(t\right)}{r}\\ m\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }y\left(t\right)=\frac{F\left(t\right)\hspace{0.17em}y\left(t\right)}{r}-g\hspace{0.17em}m\\ {x\left(t\right)}^2+{y\left(t\right)}^2=r^2\end{array}\right)$$

vars = [x(t), y(t), F(t)]

vars = $$\left(\begin{array}{ccc}x\left(t\right)& y\left(t\right)& F\left(t\right)\end{array}\right)$$

Rewrite this DAE system to a system of first-order differential algebraic equations.

[eqs, vars, newVars] = reduceDifferentialOrder(eqs, vars)

eqs = 
$$\left(\begin{array}{c}m\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }\mathrm{Dxt}\left(t\right)-\frac{F\left(t\right)\hspace{0.17em}x\left(t\right)}{r}\\ g\hspace{0.17em}m+m\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }\mathrm{Dyt}\left(t\right)-\frac{F\left(t\right)\hspace{0.17em}y\left(t\right)}{r}\\ -r^2+{x\left(t\right)}^2+{y\left(t\right)}^2\\ \mathrm{Dxt}\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }x\left(t\right)\\ \mathrm{Dyt}\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }y\left(t\right)\end{array}\right)$$

vars = 
$$\left(\begin{array}{c}x\left(t\right)\\ y\left(t\right)\\ F\left(t\right)\\ \mathrm{Dxt}\left(t\right)\\ \mathrm{Dyt}\left(t\right)\end{array}\right)$$

newVars = 
$$\left(\begin{array}{cc}\mathrm{Dxt}\left(t\right)& \frac{\partial }{\partial t}\mathrm{ }x\left(t\right)\\ \mathrm{Dyt}\left(t\right)& \frac{\partial }{\partial t}\mathrm{ }y\left(t\right)\end{array}\right)$$

Try Solving the High-Index DAE System

Before you can use a numerical MATLAB solver, such as ode15i, you must follow these steps.

To convert a DAE system to a MATLAB function handle, use daeFunction.

F = daeFunction(eqs, vars, [m, g, r])

F = function_handle with value:
    @(t,in2,in3,in4)[in3(4,:).*in4(:,1)-(in2(3,:).*in2(1,:))./in4(:,3);in3(5,:).*in4(:,1)+in4(:,1).*in4(:,2)-(in2(3,:).*in2(2,:))./in4(:,3);-in4(:,3).^2+in2(1,:).^2+in2(2,:).^2;in2(4,:)-in3(1,:);in2(5,:)-in3(2,:)]

Assign numerical values to the symbolic parameters of the system: m = 1kg, g = 9.18m/s^2, and r = 1m.

f = @(t, y, yp)  F(t, y, yp, [1, 9.81, 1])

f = function_handle with value:
    @(t,y,yp)F(t,y,yp,[1,9.81,1])

The function handle f is a suitable input for the numerical solver ode15i. The next step is to compute consistent initial conditions. Use odeset to set numerical tolerances. Then use the MATLAB decic function to compute consistent initial conditions y0, yp0 for the positions and the derivatives at time t0 = 0.

opt = odeset('RelTol', 10.0^(-4), 'AbsTol' , 10.0^(-4));
t0 = 0;
[y0,yp0] = decic(f, t0, [0.98;-0.21; zeros(3,1)], [], zeros(5,1), [], opt)

y0 = 5×1

    0.9777
   -0.2100
         0
         0
         0

yp0 = 5×1

         0
         0
         0
         0
   -9.8100

Test the initial conditions:

f(t0, y0, yp0)

ans = 5×1
10-16 ×

         0
         0
   -0.3469
         0
         0

Now you can use ode15i to try solving the system. When you call ode15i, the integration stops immediately and issues the following warnings.

Warning: Matrix is singular, close to singular or badly scaled.

Results may be inaccurate. RCOND = NaN.

Warning: Failure at t=0.000000e+00.

Unable to meet integration tolerances without reducing the step

size below the smallest value allowed (0.000000e+00) at time t.

For this example, ode15i issues these warnings multiple times. For readability, disable warnings by using warning('off','all') before calling ode15i and then enable them again.

tfinal = 0.5;
s = warning('off','all');
ode15i(f, [t0, tfinal], y0, yp0, opt);

warning(s)

Analyze and Adjust the DAE System

Check the differential index of the DAE system.

isLowIndexDAE(eqs, vars)

ans = logical
   0

This result explains why ode15i cannot solve this system. This function requires the input DAE system to be of differential index 0 or 1. Reduce the differential index by extending the model to an equivalent larger DAE system that includes some hidden algebraic constraints.

[eqs, vars, newVars, index] = reduceDAEIndex(eqs, vars)

eqs = 
$$\left(\begin{array}{c}m\hspace{0.17em}\mathrm{Dxtt}\left(t\right)-\frac{F\left(t\right)\hspace{0.17em}x\left(t\right)}{r}\\ g\hspace{0.17em}m+m\hspace{0.17em}\mathrm{Dytt}\left(t\right)-\frac{F\left(t\right)\hspace{0.17em}y\left(t\right)}{r}\\ -r^2+{x\left(t\right)}^2+{y\left(t\right)}^2\\ \mathrm{Dxt}\left(t\right)-{\mathrm{Dxt}}_1\left(t\right)\\ \mathrm{Dyt}\left(t\right)-{\mathrm{Dyt}}_1\left(t\right)\\ 2\hspace{0.17em}{\mathrm{Dxt}}_1\left(t\right)\hspace{0.17em}x\left(t\right)+2\hspace{0.17em}{\mathrm{Dyt}}_1\left(t\right)\hspace{0.17em}y\left(t\right)\\ 2\hspace{0.17em}y\left(t\right)\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }{\mathrm{Dyt}}_1\left(t\right)+2\hspace{0.17em}{{\mathrm{Dxt}}_1\left(t\right)}^2+2\hspace{0.17em}{{\mathrm{Dyt}}_1\left(t\right)}^2+2\hspace{0.17em}\mathrm{Dxt1t}\left(t\right)\hspace{0.17em}x\left(t\right)\\ \mathrm{Dxtt}\left(t\right)-\mathrm{Dxt1t}\left(t\right)\\ \mathrm{Dytt}\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }{\mathrm{Dyt}}_1\left(t\right)\\ {\mathrm{Dyt}}_1\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }y\left(t\right)\end{array}\right)$$

vars = 
$$\left(\begin{array}{c}x\left(t\right)\\ y\left(t\right)\\ F\left(t\right)\\ \mathrm{Dxt}\left(t\right)\\ \mathrm{Dyt}\left(t\right)\\ \mathrm{Dytt}\left(t\right)\\ \mathrm{Dxtt}\left(t\right)\\ {\mathrm{Dxt}}_1\left(t\right)\\ {\mathrm{Dyt}}_1\left(t\right)\\ \mathrm{Dxt1t}\left(t\right)\end{array}\right)$$

newVars = 
$$\left(\begin{array}{cc}\mathrm{Dytt}\left(t\right)& \frac{\partial }{\partial t}\mathrm{ }\mathrm{Dyt}\left(t\right)\\ \mathrm{Dxtt}\left(t\right)& \frac{\partial }{\partial t}\mathrm{ }\mathrm{Dxt}\left(t\right)\\ {\mathrm{Dxt}}_1\left(t\right)& \frac{\partial }{\partial t}\mathrm{ }x\left(t\right)\\ {\mathrm{Dyt}}_1\left(t\right)& \frac{\partial }{\partial t}\mathrm{ }y\left(t\right)\\ \mathrm{Dxt1t}\left(t\right)& \frac{{\partial }^2}{\partial t^2}\mathrm{ }x\left(t\right)\end{array}\right)$$

index = 
3

The fourth output shows that the differential index of the original model is three. Simplify the new system.

[eqs, vars, S] = reduceRedundancies(eqs, vars)

eqs = 
$$\left(\begin{array}{c}-\frac{F\left(t\right)\hspace{0.17em}x\left(t\right)-m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dxtt}\left(t\right)}{r}\\ \frac{g\hspace{0.17em}m\hspace{0.17em}r-F\left(t\right)\hspace{0.17em}y\left(t\right)+m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dytt}\left(t\right)}{r}\\ -r^2+{x\left(t\right)}^2+{y\left(t\right)}^2\\ 2\hspace{0.17em}\mathrm{Dxt}\left(t\right)\hspace{0.17em}x\left(t\right)+2\hspace{0.17em}\mathrm{Dyt}\left(t\right)\hspace{0.17em}y\left(t\right)\\ 2\hspace{0.17em}{\mathrm{Dxt}\left(t\right)}^2+2\hspace{0.17em}{\mathrm{Dyt}\left(t\right)}^2+2\hspace{0.17em}\mathrm{Dxtt}\left(t\right)\hspace{0.17em}x\left(t\right)+2\hspace{0.17em}\mathrm{Dytt}\left(t\right)\hspace{0.17em}y\left(t\right)\\ \mathrm{Dytt}\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }\mathrm{Dyt}\left(t\right)\\ \mathrm{Dyt}\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }y\left(t\right)\end{array}\right)$$

vars = 
$$\left(\begin{array}{c}x\left(t\right)\\ y\left(t\right)\\ F\left(t\right)\\ \mathrm{Dxt}\left(t\right)\\ \mathrm{Dyt}\left(t\right)\\ \mathrm{Dytt}\left(t\right)\\ \mathrm{Dxtt}\left(t\right)\end{array}\right)$$

S = struct with fields:
      solvedEquations: [3x1 sym]
    constantVariables: [0x2 sym]
    replacedVariables: [3x2 sym]
       otherEquations: [0x1 sym]

Check if the new system has a low differential index (0 or 1).

isLowIndexDAE(eqs, vars)

ans = logical
   1

Solve the Low-Index DAE System

Generate a MATLAB function handle that replaces the symbolic parameters by numerical values.

F = daeFunction(eqs, vars, [m, g, r]) 

F = function_handle with value:
    @(t,in2,in3,in4)[-(in2(3,:).*in2(1,:)-in2(7,:).*in4(:,1).*in4(:,3))./in4(:,3);(-in2(3,:).*in2(2,:)+in2(6,:).*in4(:,1).*in4(:,3)+in4(:,1).*in4(:,2).*in4(:,3))./in4(:,3);-in4(:,3).^2+in2(1,:).^2+in2(2,:).^2;in2(4,:).*in2(1,:).*2.0+in2(5,:).*in2(2,:).*2.0;in2(7,:).*in2(1,:).*2.0+in2(6,:).*in2(2,:).*2.0+in2(4,:).^2.*2.0+in2(5,:).^2.*2.0;in2(6,:)-in3(5,:);in2(5,:)-in3(2,:)]

f = @(t, y, yp)  F(t, y, yp, [1, 9.81, 1])

f = function_handle with value:
    @(t,y,yp)F(t,y,yp,[1,9.81,1])

Compute consistent initial conditions for the index reduced by the MATLAB decic function. Here, opt is the options structure that sets numerical tolerances. You already computed it using odeset.

[y0,yp0] = decic(f, t0, [0.98;-0.21; zeros(5,1)], [], zeros(7,1), [], opt)

y0 = 7×1

    0.9779
   -0.2093
   -2.0528
    0.0000
         0
   -9.3804
   -2.0074

yp0 = 7×1

         0
         0
         0
         0
   -9.3804
         0
         0

Solve the system and plot the solution.

ode15i(f, [t0, tfinal], y0, yp0, opt)