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

  • 1
    $\begingroup$ 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. $\endgroup$
    – Artes
    Mar 28, 2017 at 12:28
  • $\begingroup$ @Artes Thank you very much for the suggestion. Is it possible to see why this is the only possible solution to the third equation? $\endgroup$
    – Lukas
    Mar 28, 2017 at 13:26
  • 1
    $\begingroup$ 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. $\endgroup$
    – Artes
    Mar 28, 2017 at 14:42
  • $\begingroup$ @Artes Thanks alot for the explanation $\endgroup$
    – Lukas
    Mar 29, 2017 at 9:32

1 Answer 1


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) + 
      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}]

constraint error

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]

constraint error

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]

constraint error

  • 1
    $\begingroup$ 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 $\endgroup$
    – Lukas
    Mar 29, 2017 at 9:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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