eqns =
As per my understanding, you are trying to solve the pendulum system given in:
by using ode15s, instead of ode15i. ode15s is a specialized mass matrix solver.
This is explained in detailed way in the documentation:
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)
%Reduce the differential index of the DAEs described by eqns and vars.
[DAEs,DAEvars] = reduceDAEIndex(eqns,vars)
%Often, reduceDAEIndex introduces redundant equations and variables that can be eliminated.
% Eliminate redundant equations and variables using reduceRedundancies.
[DAEs,DAEvars] = reduceRedundancies(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)
%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)
%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 !
