MathJax TeX Test Page A more informative example is the modified double pendula problem (mod2pend.m), which is an arbitrarily made-up DAE for illustrating the block upper triangular form. \begin{align*} 0 = f_1 &= x'' + x\lambda\\[1ex] 0 = f_2 &= y'' + y\lambda + (x')^3 - G \\[1ex] 0 = f_3 &= x^2 + y^2 - L^2 \\[1ex] 0 = f_4 &= u'' + u\mu \\[1ex] 0 = f_5 &= (v''')^3 + v\mu - G \\[1ex] 0 = f_6 &= u^2 + v^2 - (L+ \alpha\lambda)^2 + \lambda'' \end{align*} Here $x, y, \lambda$ are the state variables for the first pendulum, and $u, v, \mu$ are for the second one. The first pendulum influences the second one through $\lambda$, but is not influenced by the latter. $G$ is gravity, $L>0$ is the length of both pendula, and $\alpha$ is a constant.

Visualizing structure of a DAE by showStruct()

















​​The translation into a Matlab function is mod2pend.m:


function f = mod2pend(t, z, G, L, alpha)
x = z(1); y = z(2); lam = z(3);
u = z(4); v = z(5);  mu = z(6);

f(1) = Dif(x,2) + x*lam;
f(2) = Dif(y,2) + y*lam + Dif(x)^3 - G;
f(3) = x^2 + y^2 - L^2;
f(4) = Dif(u,2) + u*mu;
f(5) = Dif(v,3)^3 + v*mu - G;
f(6) = u^2+v^2 - (L+alpha*lam)^2 + Dif(lam,2);


The script sa_mod2pend.m below displays the original structure of this DAE and its block upper triangular form.


G = 9.81; L = 10; alpha = 0.1; n = 6;

sadata = daeSA(@mod2pend, n, G, L, alpha);

figure, showStruct(sadata)

figure, showStruct(sadata, 'disptype', 'fineblocks')

© Gary Guangning Tan, 2015