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

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 and h==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) /. sol1, (h[F, delp] - 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. *)

• What is a difference with previous post you deleted on? Apr 12 at 3:59
• Number of equations is 10, number of variables including F, delp is 12, number of boundary conditions is 7. We can't solve this problem with using NDSolve and FindRoot only. Maybe we can solve it with NMinimize. Apr 12 at 6:15
• It is the same set of equations with a little bit different approach @Alex Trounev. I guess other equations are algebraic and can be solved by using another method as you mentioned but do not exactly how. My mathematica knowledge is elementary. Apr 12 at 11:05
• Solution with NMinimize is very far from True. It looks like there is no solution for this problem. Apr 12 at 11:24
• Thanks @Alex Trounev. Apr 12 at 19:52

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[%[], 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[]}, Simplify[eq[[2 ;;]] /. r'[t] -> -h[t] u[t]]];
ReplacePart[%, 6 -> Simplify[%[]/h[t]]];
eq1 = ReplacePart[%, 4 -> Simplify[%[] r[t] + %[] 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[] /.
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, h == 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 by v == 1.

bc1 = {r[1/5] == 11/5, u[1/5] == 0, v[1/5] == 19, r == 1, v == 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, v}, {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 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.

• @Dibbo123 No, you are not short on boundary conditions. However, I believe that the particular choice of values in your question does not admit a solution. For instance, a smaller value of Br` might work better. Apr 16 at 4:28