Example in the DAE solver page breaks down

The file solve-differential-algebraic-equations.html distributed with Matlab 2019b describes a step-by-step procedure for solving the equations of motion of a simple pendulum formulated as a constrained system of differential-algebraic equations. You don't need to type in the code yourself --- the Open Live Script button on that page opens a worksheet with code that you can execute.

Ultimately it solves the system of DAEs and plots the solution on the time interval [0.0 0.5].

That's about as far as the solution can be extended. Changing the 0.5 to 0.6 we get an error message saying that the solution cannot be extended beyond 0.5063. The help page makes no mention of this. I tried playing around with the tolerances but did not get anywhere.

This is puzzling since the equations model the periodic oscillations of a pundulum the solutions to which are defined for all times. The DAE solver does not even reach a small fraction of a period before it fails. Can someone shed light on what is going on?

[Steven Lord]

What is the full and exact text of the error message you received when you changed the time span in the example? Show all of the text displayed in red (and if there are any warning messages shown in orange, show that text too.)

[Rouben Rostamian]

It's not an error message -- it's a warning in orange. It says:

Warning: Failure at t=5.063260e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.798831e-15) at time t.

I am new to Matlab and therefore don't know what to do at this point.

[DC Fang]

I have exactly the same problem regarding changing the time intervals. My guess is that the constraint drifted away from x^2+y^2=r^2 eventually causing the DAE solver to fail the constraint test. Mathematica has demonstrated the same pendulum DAE problem https://reference.wolfram.com/language/tutorial/NDSolveDAE.html#1966826857 and its NDsolve function uses the IDA solver from SUNDIALS which eventually uses its covde BDF solver, and the results seem much more promising as the constraint drift is much slower.

Answer 1

[moebius]

I'm trying to solve the problem by dividing the time interval into smaller ones and using decic at the begin of each time interval. Some problems remain at the point of maximum tension T

close all

clear all

syms x(t) y(t) T(t) m r g dotx doty dotT X Y DXT Dxt(t) DYT Dyt(t) dotDxt dotDyt TT

mass=1;

radius=1;

gravity=9.81;

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

vars = [x(t); y(t); T(t)];

origVars = length(vars);

%%%% Find the equilibrium conditions %%%%%%

EqEqus = subs(eqns,{diff(x(t), 2),diff(y(t), 2),x(t),y(t),T(t),m,r,g},{0,0,X,Y,TT,mass,radius,gravity});

S=solve(EqEqus,[X,Y,TT]);

%%%%%%%%%%%%%%%%%%%%%%%

INC_M = incidenceMatrix(eqns,vars)

%M = incidenceMatrix(eqns,vars)

%%%% Reduces the order of differential equations to the lowest %%%%

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

if(isLowIndexDAE(eqns,vars)==0)

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

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

end

if(isLowIndexDAE(DAEs,DAEvars)==0)

print('problema')

end

ODEs = reduceDAEToODE(eqns, vars)

[M,f] = massMatrixForm(ODEs,vars)

pDAEs = symvar(ODEs);

pDAEvars = symvar(vars);

extraParams = setdiff(pDAEs,pDAEvars)

M = odeFunction(M, vars, g, m, r);

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

g = 9.81;

m = 1;

r = 1;

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

M = @(t, Y) M(t, Y, g, m, r);

y0est = [r*sin(pi/6); -r*cos(pi/6); m*g; 0; 0];

yp0est = zeros(5,1);

tsolold=0

for i=1:10

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

implicitODE = @(t,y,yp) M(t,y)*yp - F(t,y);

[y0, yp0] = decic(implicitODE, 0, y0est, [0 0 0 0 0], yp0est, [0 0 0 0 0], opt)

opt = odeset(opt, 'InitialSlope', yp0,'AbsTol' , 10.0^(-9),'RelTol', 0.0000000001);

[tSol,ySol] = ode23t(F, [tsolold,tsolold+0.01], y0, opt);

figure(1)

plot(tSol,ySol(:,1:origVars),'-o');

hold on

figure(2)

plot(tSol,ySol(:,1).^2+ySol(:,2).^2-r^2);

hold on

y0est=[ySol(end,1),ySol(end,2),ySol(end,3),ySol(end,4),ySol(end,5)]';

yp0est=[ySol(end,4),ySol(end,5),0,0,0]';

tsolold=tSol(end)

end