solve the pendulum problem using ode15s

Hi guys,

The pendulum problem at

https://ww2.mathworks.cn/help/symbolic/solve-differential-algebraic-equations.html?s_tid=srchtitle_site_search_1_Solve%20Differential%20Algebraic%20Equations%20%2528DAEs%2529

has been solved using ode15i.

Now I am curious that if the DAEs can be resolved with ode15s. If so, how to format the codes?

Cheers

Tony

Answer 1

[Shivam Gothi]

Hello @Tony Cheng,

As per my understanding, you are trying to solve the pendulum system given in:

https://in.mathworks.com/help/symbolic/solve-differential-algebraic-equations.html?s_tid=srchtitle_site_search_1_Solve%20Differential%20Algebraic%20Equations%20%2528DAEs%2529&requestedDomain=in

by using ode15s, instead of ode15i. ode15s is a specialized mass matrix solver.

This is explained in detailed way in the documentation:

https://in.mathworks.com/help/symbolic/solve-daes-using-mass-matrix-solvers.html

I am combining the codes given in both the above mentioned article and attaching it for your reference. It solves the pendulum system by using mass-matrix solvers.

%Specify independent variables and state variables by using syms.

syms x(t) y(t) T(t) m r g

%Specify equations by using the == operator.
eqn1 = m*diff(x(t), 2) == T(t)/r*x(t);
eqn2 = m*diff(y(t), 2) == T(t)/r*y(t) - m*g;
eqn3 = x(t)^2 + y(t)^2 == r^2;
eqns = [eqn1 eqn2 eqn3];

%Place the state variables in a column vector. Store the number of original variables for reference.
vars = [x(t); y(t); T(t)];
origVars = length(vars);

%The differential order of a DAE system is the highest differential order of its equations. 
% To solve DAEs using MATLAB, the differential order must be reduced to 1. 
% Here, the first and second equations have second-order derivatives of x(t) and y(t). 
% Thus, the differential order is 2.
%Reduce the system to a first-order system by using reduceDifferentialOrder. 
% The reduceDifferentialOrder function substitutes derivatives with new variables, 
% such as Dxt(t) and Dyt(t). The right side of the expressions in eqns is 0.

[eqns,vars] = reduceDifferentialOrder(eqns,vars)

eqns = 

vars = 

%Reduce the differential index of the DAEs described by eqns and vars.

[DAEs,DAEvars] = reduceDAEIndex(eqns,vars)

DAEs = 

DAEvars = 

%Often, reduceDAEIndex introduces redundant equations and variables that can be eliminated. 
% Eliminate redundant equations and variables using reduceRedundancies.

[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars)

DAEs = 

DAEvars = 

%%%%%%%%%%%%%%%%%%%%%%%%%%%% ODE15s %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%From the output of reduceDAEIndex, you have a vector of equations DAEs and a vector of variables DAEvars. 
% To use ode15s or ode23t, 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 where M(t,y(t))y’(t) = f(t,y(t))

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

M = 

f = 

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

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

extraParams = 

%The mass matrix M does not have these extra parameters. Therefore, convert M directly 
% to a function handle by using odeFunction.

M = odeFunction(M, DAEvars);

%Convert f to a function handle. Specify the extra parameters as additional inputs to odeFunction.

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

%Find Initial Conditions
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 matrix M and initial guesses yp0est, and 
% specifies numerical tolerances for the numerical search.
opt = odeset('Mass', M, 'InitialSlope', yp0est);

%Find consistent initial values for the variables and their derivatives by using the 
% MATLAB decic function. The first argument of decic must be a function handle describing 
% the DAE as f(t,y,yp) = f(t,y,y') = 0. In terms of M and F, this means 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);

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

%Solve DAE System
%Solve the system integrating over the time span 0 ≤ t ≤ 0.5. Add the grid lines and the legend 
% to the plot. The code uses ode15s. Instead, you can use ode23s by replacing ode15s with ode23s.

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

for k = 1:origVars-1
  S{k} = char(DAEvars(k));
end

legend(S, 'Location', 'Best');
xlabel('time(s)');
grid on;

figure;
plot(tSol,ySol(:,3)); %plots the T(t) on different graph
title("T(t)");
xlabel("time(s)")

I hope this helps !

[Tony Cheng]

Hi Shivam,

Many thanks for your nice and great help! I am trying to learn your idea in the codes and compare them with mine.

However, in my codes, the variables are actually the first five variables as in your DAEvars, and then the size of the mass matrix is 5 * 5. I am not sure whether my codes are right or not.

I will compute the numerical results and then compare them with yours. Several days later I will post my codes after I check them.