equation solving - NDSolve solution violates initial conditions

I have the following code:

solution = NDSolve[{5269.333333333333` Cos[a[t]] + 1.` Cos[a[t]] l[t] + 
 83.33333333333333` Cos[a[t] - c[t]] Derivative[1][c][t]^2 + 
 172.66666666666666` Derivative[2][a][t] + 
 8.` Cos[a[t] + b[t]] Derivative[2][b][t] == 
8.` Sin[a[t] + b[t]] Derivative[1][b][t]^2 + 
 83.33333333333333` Sin[a[t] - c[t]] Derivative[2][c][t], 
 Cos[b[t]] l[t] + 6 Cos[a[t] + b[t]] Derivative[2][a][t] + 
 8 Derivative[2][b][t] == 

64 Cos[b[t]] + 6 Sin[a[t] + b[t]] Derivative[1][a][t]^2, 
 1.` Sin[c[t]] + 0.015625` Derivative[2][c][t] == 
0.03125` Cos[a[t] - c[t]] Derivative[1][a][t]^2 + 
 0.03125` Sin[a[t] - c[t]] Derivative[2][a][t], 
6 Sin[a[t]] + 8 Sin[b[t]] == 3 Sqrt[2], 
4 a[0] == \[Pi], b[0] == 0, c[0] == 0, 
Derivative[1][a][0] == 0, 
Derivative[1][b][0] == 0, 
Derivative[1][c][0] == 0},
{a[t], b[t], c[t], l[t]},

{t, 0., 0.25}, Method -> {"IndexReduction" -> Automatic}];

asol[t_] = a[t] /. Flatten[solution];
Print["a[0]=", asol[0] , "= and a'[t]=", Derivative[1][asol][0]]

Note that I have a'[t] = Derivative[1][a][0] == 0 among the initial conditions. Yet, the output of this cell is

a[0]=0.785398= and a'[t]=6.06109

a'[t] != 0! I tried restarting Mathematica and pasting this into a new notebook, same thing. When I plot a[t], it indeed trends up instead of starting with a slope of 0. I suspect the odds of me discovering a bug in NDSolve the first time I use it are about 0 (or 0.`5) so I suspect I am not using it right.

What am I doing wrong here? Why is Mathematica giving me a solution that is NOT a solution? Any pointer appreciated.

Answer

Use the method option

Method -> {"IndexReduction" -> {Automatic, "ConstraintMethod" -> "Projection"}}

This forces the equations to be incorporated as constraints. See tutorial/NDSolveDAE#128085219. Depending on the version, you might need to us Rationalize to make the coefficients exact to avoid 1/0 errors. (In general, I avoid machine precision coefficients when doing algebra, especially in a case like this where there's numerical inconsistency. Full code below.)

With this setting I get the following:

Print["a[0]=", asol[0], "= and a'[t]=", Derivative[1][asol][0]]

a[0]=0.785398= and a'[t]=-2.77556*10^-17

---

Update: Code dump

solution = 
  NDSolve[Rationalize@{5269.333333333333` Cos[a[t]] + 
       1.` Cos[a[t]] l[t] + 
       83.33333333333333` Cos[a[t] - c[t]] Derivative[1][c][t]^2 + 
       172.66666666666666` Derivative[2][a][t] + 
       8.` Cos[a[t] + b[t]] Derivative[2][b][t] == 
      8.` Sin[a[t] + b[t]] Derivative[1][b][t]^2 + 

       83.33333333333333` Sin[a[t] - c[t]] Derivative[2][c][t], 
     Cos[b[t]] l[t] + 6 Cos[a[t] + b[t]] Derivative[2][a][t] + 
       8 Derivative[2][b][t] == 
      64 Cos[b[t]] + 6 Sin[a[t] + b[t]] Derivative[1][a][t]^2, 
     1.` Sin[c[t]] + 0.015625` Derivative[2][c][t] == 
      0.03125` Cos[a[t] - c[t]] Derivative[1][a][t]^2 + 
       0.03125` Sin[a[t] - c[t]] Derivative[2][a][t], 
     6 Sin[a[t]] + 8 Sin[b[t]] == 3 Sqrt[2], 4 a[0] == \[Pi], 
     b[0] == 0, c[0] == 0, Derivative[1][a][0] == 0, 
     Derivative[1][b][0] == 0, Derivative[1][c][0] == 0}, {a[t], b[t],

     c[t], l[t]}, {t, 0., 0.25}, 
   Method -> {"IndexReduction" -> {Automatic, 
       "ConstraintMethod" -> "Projection"}}];

asol[t_] = a[t] /. Flatten[solution];
Print["a[0]=", asol[0], "= and a'[t]=", Derivative[1][asol][0]]