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

Question

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

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