Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question

1 Answer 1

up vote 5 down vote accepted

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]]
share|improve this answer
    
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 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 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 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 at 18:10

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.