Step 1. Convert DAEs to Function Handles

From the output ofreduceDAEIndex, you have a vector of equationsDAEsand a vector of variablesDAEvars. To useode15sorode23t, 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 whereM(t,y(t))y’(t) =f(t,y(t)).

Find these function handles by usingmassMatrixForm质量矩阵Mand the right sidesF.

[M,f] = massMatrixForm(DAEs,DAEvars)

M = [ 0, 0, 0, 0, 0, 0, 0] [ 0, 0, 0, 0, 0, 0, 0] [ 0, 0, 0, 0, 0, 0, 0] [ 0, 0, 0, 0, 0, 0, 0] [ 0, 0, 0, 0, 0, 0, 0] [ 0, 0, 0, 0, -1, 0, 0] [ 0, -1, 0, 0, 0, 0, 0] f = (T(t)*x(t) - m*r*Dxtt(t))/r -(g*m*r - T(t)*y(t) + m*r*Dytt(t))/r r^2 - y(t)^2 - x(t)^2 - 2*Dxt(t)*x(t) - 2*Dyt(t)*y(t) - 2*Dxtt(t)*x(t) - 2*Dytt(t)*y(t) - 2*Dxt(t)^2 - 2*Dyt(t)^2 -Dytt(t) -Dyt(t)

The equations inDAEscan contain symbolic parameters that are not specified in the vector of variablesDAEvars. Find these parameters by usingsetdiffon the output ofsymvarfromDAEsandDAEvars.

pDAEs = symvar(DAEs); pDAEvars = symvar(DAEvars); extraParams = setdiff(pDAEs, pDAEvars)

extraParams = [ g, m, r]

The mass matrixMdoes not have these extra parameters. Therefore, convertMdirectly to a function handle by usingodeFunction.

M = odeFunction(M, DAEvars);

Convertfto a function handle. Specify the extra parameters as additional inputs toodeFunction.

f = odeFunction(f, DAEvars, g, m, r);

The rest of the workflow is purely numerical. Set parameter values and create the function handle.

g = 9.81; m = 1; r = 1; F = @(t, Y) f(t, Y, g, m, r);

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 MATLABdecicfunction. Thedecicaccepts guesses (which might not satisfy the equations) for the initial conditions, and tries to find satisfactory initial conditions using those guesses.deciccan fail, in which case you must manually supply consistent initial values for your problem.

First, check the variables inDAEvars.

DAEvars

DAEvars = x(t) y(t) T(t) Dxt(t) Dyt(t) Dytt(t) Dxtt(t)

Here,Dxt(t)is the first derivative ofx(t),Dytt(t)is the second derivative ofy(t), and so on. There are 7 variables in a7-by-1vector. Thus, guesses for initial values of variables and their derivatives must also be7-by-1vectors.

Assume the initial angular displacement of the pendulum is 30° orpi/6, and the origin of the coordinates is at the suspension point of the pendulum. Given that we used a radiusrof1, the initial horizontal positionx(t)isr*sin(pi/6). The initial vertical positiony(t)is-r*cos(pi/6). Specify these initial values of the variables in the vectory0est.

Arbitrarily set the initial values of the remaining variables and their derivatives to0. These are not good guesses. However, they suffice for our problem. In your problem, ifdecicerrors, then provide better guesses and refer to thedecicpage.

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

Create an option set that contains the mass matrixMand initial guessesyp0est, and specifies numerical tolerances for the numerical search.

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

Find consistent initial values for the variables and their derivatives by using the MATLABdecicfunction. The first argument ofdecicmust be a function handle describing the DAE asf(t,y,yp) = f(t,y,y') = 0. In terms ofMandF,这意味着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 = 0.4771 -0.8788 -8.6214 0 0.0000 -2.2333 -4.1135 yp0 = 0 0.0000 0 0 -2.2333 0 0

Now create an option set that contains the mass matrixMof the system and the vectoryp0of 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 span0≤t≤0.5. Add the grid lines and the legend to the plot. The code usesode15s. Instead, you can useode23sby replacingode15swithode23s.

[tSol,ySol] = ode15s(F, [0, 0.5], y0, opt); plot(tSol,ySol(:,1:origVars),'-o')fork = 1:origVars S{k} = char(DAEvars(k));endlegend(S,'Location','Best') gridon

Solve the system for different parameter values by setting the new value and regenerating the function handle and initial conditions.

Setrto2and regenerate the function handle and initial conditions.

r = 2; F = @(t, Y) f(t, Y, g, m, r); y0est = [r*sin(pi/6); -r*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')fork = 1:origVars S{k} = char(DAEvars(k));endlegend(S,'Location','Best') gridon