IDA vs DASSL on a special case

[mottelet]

Hello, I have an old example of a pendule with sliding pivot (constrained to f(x,y)=0) in the plane which I can solve with DASSL although I have a problem with IDA, right at startup (At t = 0 and h = 5.04637e-13, the error test failed repeatedly or with |h| = hmin).

The formulation seems to be of index 1 (can DASSL solve index 2 DAEs ?), here is the code of the residual (Scilab code, but should be quite readable), where m is the tip mass and M the mass of the pivot, x, u the speed, fx is df/dx and fy is df/dy:

function res=pendg(t,y,ydot)
    x          = y(1:3); // x(3) is angle theta
    u          = y(4:6);
    lambda     = y(7);
    xp         = ydot(1:3);
    up         = ydot(4:6);
    res         = [xp-u
                  (M+m)*up(1)+m*l*cos(x(3))*up(3)-m*l*sin(x(3))*u(3).^2+lambda*fx(x(1),x(2))+k*u(1)
                  (M+m)*up(2)+m*l*sin(x(3))*up(3)+m*l*cos(x(3))*u(3).^2+(M+m)*g+lambda*fy(x(1),x(2))+k*u(2)
                  m*l*cos(x(3))*up(1)+m*l*sin(x(3))*up(2)+m*l^2*up(3)+m*g*sin(x(3))
                  -(fx(x(1),x(2))*u(1)+fy(x(1),x(2))*u(2))];
end

Has someone an explanation why IDA could fail ?

S.

[Steven-Roberts]

Hi @mottelet. Your problem appears to be of the form

x' = u M(x, u) u' = k(x, u, lambda) 0 = g(x, u)

which is index-2. ~~IDA assumes the index is at most 1 which is likely why it is failing.~~

[mottelet]

Thanks for this confirmation. I thought that DASSL had the same limitation, is it the case ?

[Steven-Roberts]

Correction: IDA can solve DAEs with index>1, but some care is needed. You cannot use IDACalcIC at high indicies, so you should provide full initial conditions which satisfy F(t_0, y_0, y'_0) = 0. Also, you may want to use IDASetSuppressAlg to avoid algebraic errors being estimated inaccurately.

[Steven-Roberts]

Could you share the tolerances, linear solver, nonlinear solver, and other relevant configuration you are using? As far as I can tell, your problem should be solvable by IDA (and DASSL).

I put together a test to attempt to replicate the error you are seeing. It solves the slightly simpler pendulum problem

x_1'=u_1 x_2'=u_2 u_1'=-x_1 * lambda u_2'=-x_2 * lambda - 1 0=x_1 * u_1 + x_2 * u_2

which is also index-2. I did not encounter the error you have at t=0. If you can provide the fx and fy functions you use, I could try a test problem that's closer to your's.

pendulum.zip

[mottelet]

Thanks a lot, using IDASuppressAlg solved it for the sliding problem. For the above pendulum I compared with and without indication of algebric variable and suppressing it in the error test. My program below

function out=res(t,y,yd)
    x=y(1:2);
    xd=yd(1:2);
    u=y(3:4);
    ud=yd(3:4);
    lambda=y(5);
    out=[xd-u
         ud+x*lambda+[0;1]
         x'*u];
endfunction

function out=jac(t,y,yd,c)
    x=y(1:2);
    xd=yd(1:2);
    u=y(3:4);
    ud=yd(3:4);
    lambda=y(5);
    I=eye(2,2);
    out=[c*I -I zeros(2,1)
         lambda*I c*I x
         u' x' 0];
endfunction

tspan = [0,7.4162987092054876737354013887810];
y0=[1;0;0;0;0];
yd0=[0;0;0;-1;0];

[t,y,yd,info1]=ida(res,tspan(2),y0,yd0,t0=0,maxSteps=1e8,jacobian=jac)
mprintf("\nError: %g\nElsapsed time: %g\n",y(1)-1,info1.stats.eTime)

[t,y,yd,info2]=ida(res,tspan(2),y0,yd0,t0=0,yIsAlgebric=5,suppressAlg=%t,jacobian=jac)
mprintf("\nError: %g\nElsapsed time: %g\n",y(1)-1,info2.stats.eTime)

yields the output

Error: -0.00215389
Elsapsed time: 0.970291

Error: -5.76078e-05
Elsapsed time: 0.008635

it shows that adding all the algebric variable info help a lot.

[Steven-Roberts]

Glad that worked out. A couple notes about IDASuppressAlg, though:

• My understanding is that DASSL does not have a similar feature; it does not distinguish between the differential and algebraic variables for error estimation. So this does not explain why DASSL initially worked but IDA drove the first step to 0.
• It is still surprising the step size controller failed in this way when IDASuppressAlg was disabled. I would expect the error estimation to be slightly inaccurate but not too the point where h->0.

We can investigate this further. Perhaps there is a slight difference in the default configuration of DASSL and IDA that explains this.

[mottelet]

In the next days I will write a pure C implementation of the failing example to allow further investigations.

[mottelet]

In your pendulum.c example (and in my Scilab code), lambda is initially 0 which is correct considering the initial position (0,1), but when I test with (sqrt(2)/2,sqrt(2)/2)), initializing lambda with zero seems fine although the correct value, which is attained at first (very small) step is 6.9367175. Is it OK to initialize multipliers to zero for a DAE of this kind (comming from a Lagrangian of some pointwise masses with constraints) ?

[Steven-Roberts]

Is it OK to initialize multipliers to zero for a DAE of this kind (comming from a Lagrangian of some pointwise masses with constraints) ?

Using inconsistent conditions is risky and can cause failures like the Newton iteration diverging. See https://doi.org/10.1016/S0098-1354(00)00655-4 for examples with DASSL and other DAE solvers. If a BDF-based solver manages to complete the first step despite the inconsistent initial conditions, it will produce a new solution that satisfies the algebraic constraint (at least to the tolerances). Then it can proceed normally with consistent initial conditions.

For the pendulum problem starting at (0,1), I tried perturbing the initial lambda from 0 to 1e-5. This causes the error

At t = 0 and h = 7.26035e-12, the error test failed repeatedly or with |h| = hmin.

So the consistent initialization of lambda is important.

[mottelet]

Here is the example (a sliding pendulum) that fails with IDA: sliding_pendulum.c.zip unless IDASuppressAlg is used (the corresponding part has been commented out in the source). In that case the output is

y[0]=0.994789
y[1]=0.659592
y[2]=0.0052094
y[3]=-0.209267
y[4]=-0.448751
y[5]=0.209174
y[6]=-2.04045

The curve has the non-parametric equation y = x^2+cos(omega*x)/3. Note: IDA also succeeds without calling IDASuppressAlg when the Lagrange multiplier is correctly initialized (with value -1.9141723).

[mottelet]

Here is the same example with the code using dassl : sliding_pendulum_dassl.c.zip and default options. I got the following output (at t=0.05):

y[0]=0.994788
y[1]=0.659589
y[2]=0.00521073
y[3]=-0.209267
y[4]=-0.44875
y[5]=0.209173
y[6]=-2.0401

[Steven-Roberts]

@mottelet Thanks for setting up theses C examples. I'll get them running and see if I can understand the discrepancy.

[mottelet]

The dassl example has been linked against netlib's Slatec and both examples against macOS Lapack.

[Steven-Roberts]

Hi @mottelet, I was finally able to trace down the difference between IDA and DASSL to a bug in DASSL. When the first step is repeatedly rejected, it incorrectly scales some internal coefficients. This was later fixed in DASPK and is implemented correctly in IDA.

This may seem counterintuitive because DASSL completed the integration, but IDA produced an error. The real issue here is the initial y_0 and y'_0. While they are consistent in the sense that F(t_0, y_0, y'_0) = 0, I believe they are inconsistent in the sense y'_0 != y'(t_0), where y(t) is the exact trajectory starting from y_0. This manifests in an instantaneous jump in lambda from 0 to -1.9141723 which you identified. No matter how small the time step is, the jump will cause a large error estimate and rejected steps. With fully consistent initial conditions, both DASSL and IDA ought to run without errors.

[mottelet]

OK thanks. More generally speaking, would you recommend using IDA instead of DASLL suite in every situation ?

[Steven-Roberts]

IDA is the progression of DASSL and DASPK, and it should generally be preferred over the older incarnations. IDA offers more flexibility with the configuration, bug fixes, and is actively maintained.

[cswoodward]

Hi All. I can't speak fully to what is currently in DASSL and the more recent DASPK2.0 code (which are hosted here: https://cse.cs.ucsb.edu/software). IDA started as a C language rewrite of the very robust DASPK code. Recent additions and changes to IDA allow for new flexibility for solvers and data structures, and IDA is now released with support for several GPU-based programming models.

[mottelet]

OK. I was just asking because some software (e.g. DiffEq in Julia) propose DASSL. But the logic of DiffEq is to rewrite all solvers natively, hence simpler to do for DASSL which is standalone...