5
$\begingroup$

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.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
5
$\begingroup$

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]]
Improve this answer
$\endgroup$
4
  • $\begingroup$ 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 $\endgroup$ Commented Aug 23, 2014 at 13:21
  • $\begingroup$ @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)" $\endgroup$ Commented Aug 23, 2014 at 16:22
  • $\begingroup$ @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. $\endgroup$ Commented Aug 23, 2014 at 17:21
  • $\begingroup$ @ThomasMaterna Try using Rationalize (see update) to avoid 1/0 errors. I assume there's some round-off error somewhere. $\endgroup$ Commented Aug 23, 2014 at 18:10

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.