# Find numerical solution of (overdetermined) system of differential(-algebraic) equations

I am attempting to find numerical solutions to a system of differential equations I found. First off, these are the two relevant functions:

f[t_] := 0.2856366975757554*Cos[3.141592653589793*(0.3295 - 0.0002*Tan[1.560831139643837*Erf[0.2*(-15 + t)]])]^(1/4);
g[t_] := -(111/20) + 6*Sqrt[3]*Sqrt[Cos[3.141592653589793*(0.3295 - 0.0002*Tan[1.560831139643837*Erf[0.2*(-15 + t)]])]];
k = 2*Sqrt[3];


The actual system of differential equations is given by

eqs0 = { (* this is the actual system of ODEs *)
f[t]*Cos[k*h1[t]] - g[t]/k*Sin[k*h1[t]] + h2'[t] == 0,
g[t]/k*Sin[k*h2[t]] + h1'[t] == 0,
h1[t]*h2'[t] - h1'[t]*h2[t] == 0
};


where I need to find h1[t] and h2[t]. From the nature of the problem I safely know that h1[t]!=0 so that I may use the third equation to reduce the total number of equations, which will turn the problem into a system of differential-algebraic equations:

eqs1 = { (* here I used the last equation of eqs0 to reduce the number of equations *)
f[t]*Cos[k*h1[t]] - g[t]/k*Sin[k*h1[t]] + h2[t]/h1[t]*h1'[t] == 0,
g[t]/k*Sin[k*h2[t]] + h1'[t] == 0
};


The initial conditions for the problem are unfortunately not too clear, but it does definitely make sense to demand h2[0]==0. Apart from that, I actually know nothing about h2[t] - it is completely open how this function will look like. There is an initial condition for h1[0] (and similarly one for h1'[0]) based on physical thoughts, but I would not consider it too strong - so I came up with two possible initial conditions, given by

initial1 = { (* motivated by the physical problem, yet the condition for h1 is not mandatory *)
h1[0] == ArcTan[k*f[0]/g[0]]/k,
h2[0] == 0
};

initial2 = { (* another initial condition that might be more compatible with eqs1 *)
h1[0] == First@(h1[0] /. NDSolve[Join[{eqs1, h2[0] == 0}], {h1, h2}, {t, 0, 0}]),
h2[0] == 0
};


Attempting to run NDSolve on eqs0 of course leads to a complaint because the system is overdetermined, so I tried on eqs1:

sol = NDSolve[Join[{eqs1, initial1}], {h1, h2}, {t, 0, 30}];


NDSolve::pdord: Some of the functions have zero differential order, so the equations will be solved as a system of differential-algebraic equations.

NDSolve::index: The DAE solver failed at t = 0.. The solver is intended for index 1 DAE systems and structural analysis indicates that the DAE index is 2. The option Method->{"IndexReduction"->Automatic} may be used to reduce the index of the system.

So, as suggested by the second warning, I change the Method which does also yield the differential-algebraic hint but also a convergence issue:

sol = NDSolve[Join[{eqs1, initial1}], {h1, h2}, {t, 0, 30}, Method -> {"IndexReduction" -> Automatic}];


NDSolve::ndcf: Repeated convergence test failure at t == 0.; unable to continue.

If I switch to the second set of initial conditions (which I assume to be more compatible with the system), I receive a warning about the initial conditions. However, I am unsure how to specify initial values for the derivatives... If I went with initial1 I could make a guess for h1'[0] but that's all.

sol = NDSolve[Join[{eqs1, initial2}], {h1, h2}, {t, 0, 30}];


NDSolve::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended.

As I am by far no expert on solving such kind of equations in general and especially not in Mathematica: How would I need to proceed in order to find solutions for the two functions in question, h1[t] and h2[t]?

• Knowing that h1[t] h2'[t] - h1'[t] h2[t] == 0 (the third equation) it implies that h2[t] == c2 h1[t] where c2 is a constant. So your problem isn't especially hard to tackle with. Mar 28, 2017 at 12:28
• @Artes Thank you very much for the suggestion. Is it possible to see why this is the only possible solution to the third equation? Mar 28, 2017 at 13:26
• Check D[ h2[t]/h1[t], t], since the third equation is equivalent under general circumstances to D[ h2[t]/h1[t], t] == 0, this implies that h2[t]/h1[t] does not depend on t. Mar 28, 2017 at 14:42
• @Artes Thanks alot for the explanation Mar 29, 2017 at 9:32

That first set of initial conditions works well with the case where we simply drop the third DE.

f[t_] := 0.2856366975757554*
Cos[3.141592653589793*(0.3295 -
0.0002*Tan[1.560831139643837*Erf[0.2*(-15 + t)]])]^(1/4);
g[t_] := -(111/20) +
6*Sqrt[3]*
Sqrt[Cos[
3.141592653589793*(0.3295 -
0.0002*Tan[1.560831139643837*Erf[0.2*(-15 + t)]])]];
k = 2*Sqrt[3];


Here I remove that last DE and instead regard it as a constraint.

eqs1 = {f[t]*Cos[k*h1[t]] - g[t]/k*Sin[k*h1[t]] + h2'[t] == 0,
g[t]/k*Sin[k*h2[t]] + h1'[t] == 0};
constraint = h1[t]*h2'[t] - h1'[t]*h2[t] == 0;


We will solve the system and check afterwards whether the constraint happened to be enforced.

sol = NDSolve[Join[{eqs1, initial1}], {h1, h2}, {t, 0, 30}];


Let's see how the solutions appear.

GraphicsRow[{Plot[Evaluate[{h1[t]} /. sol[[1]]], {t, 0, 1}],
Plot[Evaluate[{h2[t]} /. sol[[1]]], {t, 0, 1}]}]


Now we visually check the error in that constraint.

Plot[constraint[[1]] /. sol[[1]], {t, 0, 30}]


One possible way to improve on this would be to use Method -> "Projection" to enforce the constraint as an invariant. To make this work I needed to get rid of explicit derivatives in the constraint equation.

newconstraint = constraint /. First[Solve[eqs1, {h1'[t], h2'[t]}]];
solProj =
NDSolve[Join[{eqs1, initial1}], {h1[t], h2[t]}, {t, 0, 30},
Method -> {"Projection", "Invariants" -> {newconstraint[[1]]}}]


The plot is a bit different now.

GraphicsRow[{Plot[Evaluate[{h1[t]} /. solProj[[1]]], {t, 0, 1}],
Plot[Evaluate[{h2[t]} /. solProj[[1]]], {t, 0, 1}]}]


And the error is almost an order of magnitude reduced.

Plot[newconstraint[[1]] /. solProj[[1]], {t, 0, 30}, PlotRange -> All]


Increasing AccuracyGoal beyond default will bring this down even more. Here is a variant that remained reasonably fast and gave another order of magnitude improvement on the constraint violation.

solProj =
NDSolve[Join[{eqs1, initial1}], {h1[t], h2[t]}, {t, 0, 30},
Method -> {"Projection", "Invariants" -> {newconstraint[[1]]}},
PrecisionGoal -> 15, AccuracyGoal -> 20]

Plot[newconstraint[[1]] /. solProj[[1]], {t, 0, 30},
PlotRange -> All]


• Thank you very much for this detailed answer. Especially the projection method as a tool to enforce the constraint was new and very useful for this kind of problem Mar 29, 2017 at 9:33