How to deal with differential algebraic equation (DAE) in ParametricNDSolve?

Question

I have the following system of 10 differential equations where the first two equations are algebraic for u[t] and srz[t]. We have boundary conditions r[0.2] == Br, r'[0.2] == 0.0, h[0.2] == 1.0/(Br*Dr), v[0.2] == Dr, szz[0.2] == F/(r[0.2]*h[0.2]), and other two boundary conditions r[1]==1 and h[1]==1 fix two unknown shooting parameters F and delp.

I tried the following code but got errors in my calculation.

How can I fix this? Any help would be appreciated.

ClearAll["Global`*"];

x = 1.0; Br = 2.2; Dr = 19.0;

sol1 = ParametricNDSolve[{
   r'[t] == -(u[t] h[t])/x,
   srz[t] == -x*v[t] u'[t]/h[t] + (q[t] + (v[t] r'[t] w[t])/h[t])/h[t],
   ssi[t] == 2.0*u[t]/h[t] + p[t],
   szz[t] == -(2.0*x*v[t] v'[t])/h[t] + 2.0*x*v[t] r'[t] q[t]/h[t]^2 -
      p[t],
   h[t] v[t] r'[t] + h[t] r[t] v'[t] + v[t] r[t] h'[t] == 0,
   w[t] h[t] + r'[t] q[t]*x*v[t] + (u[t] h[t]^2)/r[t] - 
     v'[t]*x*v[t] h[t] == 0,
   h[t]*(2.0*w[t]/h[t] - p[t]) + r'[t]*x*v[t] srz[t] == 0, 
   h[t] srz[t] + r'[t]*x*v[t] szz[t] == 0,
   h'[t] szz[t] + h[t] szz'[t] - srz[t]/(r[t]*x*v[t]) == r'[t]*delp,
   h'[t] srz[t] + h[t] srz'[t] - 
     h[t]^2*(2.0*w[t]/h[t] - p[t] - ssi[t])/(x*v[t] r[t]) == delp,
   r[0.2] == Br, r'[0.2] == 0.0, h[0.2] == 1.0/(Br*Dr), v[0.2] == Dr, 
   szz[0.2] == F/(r[0.2]*h[0.2])
   }, {h, r, v, r', szz}, {t, 0.2, 1}, {F, delp}, 
  Method -> {"EquationSimplification" -> {Automatic, 
      "SimplifySystem" -> True}}, 
  Method -> {"IndexReduction" -> Automatic}]

FindRoot[{(r[F, delp][1] - 1) /. sol1, (h[F, delp][1] - 1) /. 
   sol1}, {{F, 1}, {delp, 1}}]

ParametricNDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions. *)

Answer 1

As written in the question, the system of equations consists of nine ODEs, one algebraic equation, two parameters, and seven boundary conditions. So, it is not surprising that Mathematica cannot solve it. However, a series of transformations can reduce the system to four ODEs, one parameter, and five boundary conditions, which Mathematica can solve with appropriate boundary conditions.

For convenience, cast the ten equations into the form lhs == 0 and work only with the left-hand sides.

x = 1; Br = 11/5; Dr = 19;
Subtract @@@ {r'[t] == -(u[t] h[t])/x, 
   srz[t] == -x*v[t] u'[t]/h[t] + (q[t] + (v[t] r'[t] w[t])/h[t])/h[t],
   ssi[t] == 2*u[t]/h[t] + p[t], 
   szz[t] == -(2*x*v[t] v'[t])/h[t] + 2*x*v[t] r'[t] q[t]/h[t]^2 - p[t], 
   h[t] v[t] r'[t] + h[t] r[t] v'[t] + v[t] r[t] h'[t] == 0, 
   w[t] h[t] + r'[t] q[t]*x*v[t] + (u[t] h[t]^2)/r[t] - v'[t]*x*v[t] h[t] == 0, 
   h[t]*(2*w[t]/h[t] - p[t]) + r'[t]*x*v[t] srz[t] == 0, 
   h[t] srz[t] + r'[t]*x*v[t] szz[t] == 0, 
   h'[t] szz[t] + h[t] szz'[t] - srz[t]/(r[t]*x*v[t]) == r'[t]*delp, 
   h'[t] srz[t] + h[t] srz'[t] - h[t]^2*(2*w[t]/h[t] - p[t] - ssi[t])/
       (x*v[t] r[t]) == delp};

First, Integrate the fifth equation to yield r[t] h[t] v[t] - c, where c is the constant of integration. Then, rescale {ssi[t], srz[t], szz[t], p[t]} by h[t].

ReplacePart[%, 5 -> Integrate[%[[5]], t] - c];
eq = Simplify[% /. {ssi -> Function[t, ssih[t]/h[t]], 
    srz -> Function[t, srzh[t]/h[t]], 
    szz -> Function[t, szzh[t]/h[t]], p[t] -> ph[t]/h[t]}];

The net effect of these transformations is to eliminate h'[t]. Then, transform the equations such that r'[t] and v'[t] appear in only one equation each, and simplify the results.

Numerator /@ Together /@ 
Join[{eq[[1]]}, Simplify[eq[[2 ;;]] /. r'[t] -> -h[t] u[t]]];
ReplacePart[%, 6 -> Simplify[%[[6]]/h[t]]];
eq1 = ReplacePart[%, 4 -> Simplify[%[[4]] r[t] + %[[6]] 2]]

(* {h[t] u[t] + r'[t], 
    -q[t] + srzh[t] + u[t] v[t] w[t] + v[t] u'[t], 
    -ph[t] + ssih[t] - 2 u[t], 
    ph[t] r[t] + 2 h[t] u[t] + r[t] (szzh[t] + 2 w[t]), 
    -c + h[t] r[t] v[t], 
    h[t] u[t] + r[t] (-q[t] u[t] v[t] + w[t] - v[t] v'[t]), 
    -ph[t] - srzh[t] u[t] v[t] + 2 w[t], 
    srzh[t] - szzh[t] u[t] v[t], 
    -srzh[t] + delp h[t]^2 r[t] u[t] v[t] + h[t] r[t] v[t] szzh'[t], 
    h[t] ph[t] + h[t] ssih[t] - delp r[t] v[t] - 2 h[t] w[t] + r[t] v[t] srzh'[t]} *)

The algebraic equations, now five in number, are isolated by

eq1nod = Cases[eq1, z_ /; ! MemberQ[z, _'[t], 4]];

and used to eliminate five variables from eq1.

seq1nod = Flatten@Solve[Thread[eq1nod == 0], {ph[t], w[t], h[t], ssih[t], srzh[t]}];
eq2 = Numerator /@ Together /@ 
DeleteCases[Simplify[eq1 /. Join[
     seq1nod[[;; 4]], {srzh -> Function[t, szzh[t] u[t] v[t]]}]], 0]

Finally, the remaining algebraic variable, q[t], is isolated and eliminated.

Simplify[eq2[[5]] /. 
   Flatten@Solve[Thread[eq2[[2 ;; 4]] == 0], {u'[t], v'[t], szzh'[t]}]];
Simplify@Flatten@Solve[% == 0, q[t]]
eq3 = Collect[Numerator /@ Together /@ Simplify[eq2[[;; 4]] /. %], _'[t], Simplify]

(* {c u[t] + r[t] v[t] r'[t], 
    -2 c^2 r[t]^2 szzh[t] u[t]^4 v[t]^3 + 2 c r[t]^4 szzh[t]^2 
    u[t]^3 v[t]^4 + c r[t]^4 szzh[t]^2 u[t]^5 v[t]^6 + 
    2 c r[t]^2 v[t] (c szzh[t] + 2 delp r[t]^2 v[t]^2) + 
    c u[t] (4 c^2 - 8 c r[t]^2 v[t] + r[t]^4 szzh[t]^2 v[t]^2) - 
    2 r[t]^2 u[t]^2 v[t] (-2 c^2 delp r[t] v[t]^2 + 
    szzh[t] (c^2 - 3 c^2 v[t]^2 + 2 r[t]^2 v[t]^4)) + 
    4 c r[t]^4 szzh[t] (-1 + u[t]^2) v[t]^4 u'[t], 
    c r[t]^4 szzh[t]^2 v[t]^2 + 4 c delp r[t]^4 u[t] v[t]^3 + 
    c r[t]^4 szzh[t]^2 u[t]^4 v[t]^6 + 4 r[t]^2 u[t]^3 v[t]^3 
    (c^2 delp r[t] + c^2 szzh[t] - r[t]^2 szzh[t] v[t]^2) + 
    2 c u[t]^2 (2 c^2 - 4 c r[t]^2 v[t] + r[t]^4 szzh[t]^2 v[t]^4) - 
    4 c r[t]^4 szzh[t] (-1 + u[t]^2) v[t]^3 v'[t], 
    u[t] (c^2 delp - r[t] szzh[t] v[t]^2) + c r[t] v[t] szzh'[t]} *)

Turn now to the boundary conditions given in the question.

bc = {r[1/5] == Br, r'[1/5] == 0, h[1/5] == 1/(Br*Dr), v[1/5] == Dr, 
    szz[1/5] == F/(r[1/5]*h[1/5]), r[1] == 1, h[1] == 1}

The first thing to notice is that szz[1/5] == F/(r[1/5]*h[1/5]) actually is not a boundary condition at all but rather the definition of F and so should be deleted for now. Next, r'[1/5] == 0 can, and should be, replaced by u[1/5] == 0. The product r[1/5] h[1/5] v[1/5] equals 1, which determines c -> 1, after which h[1/5] == 1/(Br*Dr) has been used and no longer is needed. Finally, the result r[t] h[t] v[t] == 1 permits replacing h[1] == 1 by v[1] == 1.

bc1 = {r[1/5] == 11/5, u[1/5] == 0, v[1/5] == 19, r[1] == 1, v[1] == 1};

A standard way of solving such systems involves replacing the parameter delp by a pseudo-variable, dp[t], which of course satisfies dp'[t] == 0.

Thread[Join[eq3 /. {c -> 1, delp -> dp[t]}, {dp'[t]}] == 0];
s = NDSolve[{%, bc1}, {r[t], szzh[t], u[t], v[t]}, {t, 1/5, 1}, 
    Method -> {"Shooting", "StartingInitialConditions" -> {szzh[1/5] == 
    0.1,dp[1/5] == 1}}]

for instance. Unfortunately, with this and every other "StartingInitialConditions" I have tried, it does not converge. NMinimize, however, does converge, but not to a useful answer.

{r[1/5] == 11/5, u[1/5] == 0, v[1/5] == 19, szzh[1/5] == sz0, dp[1/5] == dp0};
Thread[Join[eq3 /. {c -> 1, delp -> dp[t]}, {dp'[t]}] == 0];
pst = ParametricNDSolveValue[{%, %%}, Norm@{r[1], v[1]}, {t, 1/5, 1}, {sz0, dp0}];
NMinimize[pst[i1, i2], {i1, i2}, Method -> "NelderMead"]

(* {2.19824, {i1 -> 766.257, i2 -> 411.22}} *)

Here is the corresponding plot (r is blue, v is orange, and u is red).

The problem is that r'[t] == -u[t]/(r[t] v[t]), and u[t]/v[t] starts small and stays that way. So, r[1] is approximately 11/5 for every calculation I have tried. Perhaps, different boundary conditions would produce better results. It is, however, important to remember that there is no guarantee that some nonlinear boundary value problems even have solutions.