# 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.

-

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

-
Did you change anything else? When I replace my Method ->... with yours, I get Power::infy: Infinite expression 1/0.^2 encountered. >> and other errors – Thomas Materna Aug 23 '14 at 13:21
@Thomas Materna - What version are you using? Michael E2's solution works on my \$Version "10.0 for Mac OS X x86 (64-bit) (June 29, 2014)" – Bob Hanlon Aug 23 '14 at 16:22
@Bob Hanlon - I am using 9.0.1.0 on Windows (64-bit). It does seem to work when I try "ConstraintMethod" -> None. I just wanted to make sure Michael E2's solution worked before accepting the answer. Thanks to you both for the help. – Thomas Materna Aug 23 '14 at 17:21
@ThomasMaterna Try using Rationalize (see update) to avoid 1/0` errors. I assume there's some round-off error somewhere. – Michael E2 Aug 23 '14 at 18:10