# NDSolve returns answer not satisfying initial conditions given to it

I have a (fairly) complicated Lagrangian that I'm trying to solve numerically in Mathematica. The first part, generating the equations of motion, works fine, but I've included them here for the sake of having reproducible code

T[a_, b_, c_,adot_, bdot_, cdot_] :=1/2 8 adot^2 + 1/2 6 bdot^2 + 1/2 0.9 cdot^2 - 5.6 adot bdot Cos[a-b] +  adot cdot 0.9 Cos[a-c]
V[a_, b_, c_] := 10 Cos[a] - 55 Cos[b] -3 Cos[c]
fConstraint[a_, c_] := - 3 Cos[a] - 9 Cos[c] + 3 / Sqrt[2];

lagrange = T[a[t], b[t], c[t],adot[t], bdot[t], cdot[t]]  -  V[a[t], b[t], c[t]];
eoma = D[D[lagrange, adot[t]]/. lagrangeRules ,t] - D[lagrange, a[t]] /. lagrangeRules //Simplify;
eomb = D[D[lagrange, bdot[t]]/. lagrangeRules ,t] - D[lagrange, b[t]] /. lagrangeRules //Simplify;
eomc = D[D[lagrange, cdot[t]]/. lagrangeRules ,t] - D[lagrange, c[t]] /. lagrangeRules //Simplify;

system = {
eoma == k[t]D[fConstraint[a[t], c[t]], a[t]],
eomb == 0,
eomc == k[t]D[fConstraint[a[t], c[t]], c[t]],
fConstraint[a[t], c[t]]==0,
a[0]==\[Pi]/4, a'[0]==0,
b[0]==0, b'[0]==0,
c[0]==-\[Pi]/2, c'[0]==0} //Simplify


is a system of non-linear, coupled, 2nd degree, ordinary differential equations. The functions a[t], b[t], and c[t] are the ones I want to solve, while k[t] is a Lagrange multiplier used to make the system follow the constraint.

To solve this system I use

sola = NDSolveValue[system, a, {t, 0, 1}, Method->{IndexReduction->Automatic}];
solb = NDSolveValue[system, b, {t, 0, 1}, Method->{IndexReduction->Automatic}];
solc = NDSolveValue[system, c, {t, 0, 1}, Method->{IndexReduction->Automatic}];


The IndexReduction->Automatic option is required, without it Mathematica tells me that this is a DAE of degree 3 and that the index has to be reduced. I believe that this is something to do with the k[t] appearing in the system of equations only in that form, while the other functions also appear as first and second derivatives. Regardless, once this method is specified, Mathematica runs without throwing any warnings or errors.

If I check the solutions at t=0, I expect to recover the initial conditions specified in the last 3 lines of system

{0.785398,0.,-1.5708,0.,0.,0.}


{sola[0], solb[0], solc[0], sola'[0], solb'[0], solc'[0]} //Chop
out = {0.785398,0,-1.5708,1.12123,0,0.264276}


Which matches the initial positions but is completely wrong for the initial speeds!

Is this a bug? Am I asking too much of NDSolve here? Is there a work around?

Interestingly, apart from not matching the initial conditions, the solution appears correct. It has all the properties I'd expect it to, and even conserves the value of the associated Hamiltonian, so I know it's not complete garbage. If I plot a[t], b[t], c[t], and k[t] they also look continuous and smooth, so I don't think there's a problem with the system being ill conditioned or spiking off to infinity anywhere. This is also true if I solve them over [-1, 1] rather than just [0, 1].

My hunch is that IndexReduction is doing something bad here. If I simplify the problem to remove the constraint (set k==0 and fConstraint == 0) then everything works fine. But I don't fully understand what it does, why it's needed, or how I could get around it.

• k[0] is not determined by the initial conditions, which must be causing trouble. Not sure why there are no errors or warnings, though. Try differentiating the constraint: D[fConstraint[a[t], c[t]] == 0, t, t] (need the constraint to be second-order so that the equations determine a''[0], b''[0], c''[0] and k[0], or at least that's how it seems to me). Commented Mar 16 at 13:10
• System can be solved step by step using constraint and 3rd equation to express c[t], k[t]. See my answer. Commented Mar 16 at 17:05
• @Goofy thanks, for the suggestion. As an alternative to @Alex solution I've found that solving for k[t] in the system in addition to differentiating the constraint allows IndexReduction to be removed and then NDSolve works. As code it's system = { eoma dConstraintdc ==eomc dConstraintda, eomb == 0, D[constraint,t] == 0, a[0]==\[Pi]/4, a'[0]==0, b[0]==0, b'[0]==0, c[0]==-\[Pi]/2, c'[0]==0} //Simplify where constraint= - 3 Cos[a[t]] - 9 Cos[c[t]] + 3 / Sqrt[2]; dConstraintda = D[constraint, a[t]]; dConstraintdc = D[constraint, c[t]]; and k is found by substitution after NDSolve Commented Mar 16 at 17:34
• Although I should add that doing it that way does raise a warning: NDSolveValue: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations. But I'll take a correct answer with a warning above a silently wrong one any day! Commented Mar 16 at 17:41

We can solve this problem using Solve first as follows

T[a_, b_, c_, adot_, bdot_, cdot_] :=
1/2  8  adot^2 + 1/2  6  bdot^2 + 1/2  0.9  cdot^2 -
5.6  adot  bdot  Cos[a - b] + adot  cdot  0.9  Cos[a - c]
V[a_, b_, c_] := 10  Cos[a] - 55  Cos[b] - 3  Cos[c]
fConstraint[a_, c_] := -3  Cos[a] - 9  Cos[c] + 3/Sqrt[2];

lagrangeRules = {adot -> a', bdot -> b', cdot -> c'};

lagrange =
T[a[t], b[t], c[t], adot[t], bdot[t], cdot[t]] - V[a[t], b[t], c[t]];
eoma = D[D[lagrange, adot[t]] /. lagrangeRules, t] -
D[lagrange, a[t]] /. lagrangeRules // Simplify;
eomb = D[D[lagrange, bdot[t]] /. lagrangeRules, t] -
D[lagrange, b[t]] /. lagrangeRules // Simplify;
eomc = D[D[lagrange, cdot[t]] /. lagrangeRules, t] -
D[lagrange, c[t]] /. lagrangeRules // Simplify;

system = {eoma == k[t] D[fConstraint[a[t], c[t]], a[t]], eomb == 0,
eomc == k[t] D[fConstraint[a[t], c[t]], c[t]],
fConstraint[a[t], c[t]] == 0, a[0] == \[Pi]/4, a'[0] == 0,
b[0] == 0, b'[0] == 0, c[0] == -\[Pi]/2, c'[0] == 0} // Simplify

Solve[2  Cos[a[t]] + 6  Cos[c[t]] == Sqrt[2], c[t]]

Out[]= {{c[t] ->
ConditionalExpression[-ArcCos[1/6 (Sqrt[2] - 2 Cos[a[t]])] +
2 \[Pi] ConditionalExpression[1, \[Placeholder]],
ConditionalExpression[1, \[Placeholder]] \[Element]
Integers]}, {c[t] ->
ConditionalExpression[
ArcCos[1/6 (Sqrt[2] - 2 Cos[a[t]])] +
2 \[Pi] ConditionalExpression[1, \[Placeholder]],
ConditionalExpression[1, \[Placeholder]] \[Element] Integers]}}


From initial condition we choose c[t]= -ArcCos[1/6 (Sqrt[2] - 2 Cos[a[t]])]. Then we define derivatives

D[-ArcCos[1/6  (Sqrt[2] - 2  Cos[a[t]])], t]

Out[]= (
Sin[a[t]] Derivative[1][a][t])/(3 Sqrt[
1 - 1/36 (Sqrt[2] - 2 Cos[a[t]])^2])

D[-ArcCos[1/6  (Sqrt[2] - 2  Cos[a[t]])], t, t]

Out[]= (Cos[a[t]] Derivative[1][a][t]^2)/(
3 Sqrt[1 -
1/36 (Sqrt[2] - 2 Cos[a[t]])^2]) + ((Sqrt[2] - 2 Cos[a[t]]) Sin[
a[t]]^2 Derivative[1][a][t]^2)/(
54 (1 - 1/36 (Sqrt[2] - 2 Cos[a[t]])^2)^(3/2)) + (
Sin[a[t]] (a^\[Prime]\[Prime])[t])/(
3 Sqrt[1 - 1/36 (Sqrt[2] - 2 Cos[a[t]])^2])


and rule

rulc = {c[t] -> -ArcCos[1/6  (Sqrt[2] - 2  Cos[a[t]])],
c'[t] -> (Sin[a[t]] Derivative[1][a][t])/(
3 Sqrt[1 - 1/36 (Sqrt[2] - 2 Cos[a[t]])^2]),
c''[t] -> (Cos[a[t]] Derivative[1][a][t]^2)/(
3 Sqrt[1 -
1/36 (Sqrt[2] - 2 Cos[a[t]])^2]) + ((Sqrt[2] -
2 Cos[a[t]]) Sin[a[t]]^2 Derivative[1][a][t]^2)/(
54 (1 - 1/36 (Sqrt[2] - 2 Cos[a[t]])^2)^(3/2)) + (
Sin[a[t]] (a^\[Prime]\[Prime])[t])/(
3 Sqrt[1 - 1/36 (Sqrt[2] - 2 Cos[a[t]])^2])};


With this rule we solve

sol0 =
Solve[3 (Sin[c[t]] - 3/10  Sin[a[t] - c[t]]  Derivative[1][a][t]^2 +
3/10  Cos[a[t] - c[t]]  (a^\[Prime]\[Prime])[t] +
3/10  (c^\[Prime]\[Prime])[t]) == 9 k[t]  Sin[c[t]] /. rulc,
k[t]][[1]]

Out[]= {k[t] -> -(1/(
3 Sqrt[1 -
1/36 (Sqrt[2] - 2 Cos[a[t]])^2]))(-Sqrt[
1 - 1/36 (Sqrt[2] - 2 Cos[a[t]])^2] -
3/10 Sin[a[t] + ArcCos[1/6 (Sqrt[2] - 2 Cos[a[t]])]] Derivative[
1][a][t]^2 +
3/10 Cos[a[t] + ArcCos[1/6 (Sqrt[2] - 2 Cos[a[t]])]] (
a^\[Prime]\[Prime])[t] +
3/10 ((Cos[a[t]] Derivative[1][a][t]^2)/(
3 Sqrt[1 -
1/36 (Sqrt[2] - 2 Cos[a[t]])^2]) + ((Sqrt[2] -
2 Cos[a[t]]) Sin[a[t]]^2 Derivative[1][a][t]^2)/(
54 (1 - 1/36 (Sqrt[2] - 2 Cos[a[t]])^2)^(3/2)) + (
Sin[a[t]] (a^\[Prime]\[Prime])[t])/(
3 Sqrt[1 - 1/36 (Sqrt[2] - 2 Cos[a[t]])^2])))};


Finally we define new system and solve it as

sys = {system[[1]], system[[2]]} /. Join[sol0, rulc];sol = NDSolve[{sys, 4  a[0] == \[Pi], Derivative[1][a][0] == 0,
b[0] == 0, Derivative[1][b][0] == 0} /. Join[sol0, rulc], {a,
b}, {t, 0, 1}];


Visualization

Plot[Evaluate[{a[t], b[t], c[t] /. rulc, k[t] /. sol0} /.
sol[[1]]], {t, 0, 1}, PlotLegends -> {a, b, c, k}, Frame -> True,
FrameLabel -> Automatic]


• Thanks, that really looks like it works! I'm slightly cautious about using inverse trig functions as substitutions however. These variables are angles in a dynamical system that are expected to do full loops, is there a risk of domain co-domain problems with these functions? Also, what's up with the noise on k shortly after the start? I'm seeing the same thing with a different solution too Commented Mar 16 at 17:31
• @ScienceSnake This is not noise, this is periodic solution which turns in monotonic solution shortly, we can see it on Plot[Evaluate[{a[t], b[t], c[t] /. rulc, k[t] /. sol0} /. sol[[1]]], {t, 0, 2 10^-4}, PlotLegends -> {a, b, c, k}, Frame -> True, FrameLabel -> Automatic]. Note that all parameters has this periodic component in the begining, for example, see  Plot[Evaluate[a[t] /. sol[[1]]], {t, 0, 2 10^-4}, Frame -> True, FrameLabel -> Automatic], Commented Mar 16 at 22:16
• Also, please, check long time solution to see periodicity , for example  sol = NDSolve[{sys, 4 a[0] == [Pi], Derivative[1][a][0] == 0, b[0] == 0, Derivative[1][b][0] == 0} /. Join[sol0, rulc], {a, b}, {t, 0, 10 Pi}] Commented Mar 16 at 22:23
• The system is chaotic (in the technical sense), as it's a triple pendulum which starts with one arm resting on the ground, so it should be aperiodic in the long run. Unless I got super lucky in the random coefficients I picked for this specific example. Commented Mar 16 at 22:39
• @ScienceSnake Maybe we can compute Poincare section for this system to see chaotic dynamics? Commented Mar 16 at 23:35