Error when using NDSolve

I am trying to solve a system of differential equations as follows:

L1 = 5.0;
L2 = 0.43;
d = 0.0056;
g = 9.8;
h1 = 5.0;
h2 = 0.43;
A1 = 0.00001;
A2 = 0.001;
p0 = 100000;
v0 = 0.001;
NDSolve[{P[t] -
1*1000*10*
L2 == (1000/(2*0.0056)*y[t]^2*
L2*(0.316*(1000*y[t]*0.0056/0.00089))^(-0.25)) +
500*01.7*y[t]^2, (1000*10*L1 -
P[t]) == (1000/(2*0.0056)*z[t]^2*
L1*(0.316*(1000*z[t]*0.0056/0.00089))^(-0.25)) +
500*01.7*z[t]^2, (y[t] - z[t])*0.0001 == v'[t],
P[t] == p0*v0/v[t],P==100000,v==0.001}, {P, v, y, z}, {t, 0, 1}]

Uppon ruuning the comand I get the following error message:"NDSolve::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.". How can I make this work ?

• Is this a system of differential equations or something else? – zhk Jun 28 '18 at 1:10
• @zhk Looks like a DAE system with 1 differential equation and 3 algebraic ones. – Chris K Jun 28 '18 at 2:15

I plugged your initial conditions into your equations

{P[t] - 1*1000*10*L2 == (1000/(2*0.0056)*y[t]^2*
L2*(0.316*(1000*y[t]*0.0056/0.00089))^(-0.25)) + 500*01.7*y[t]^2,
(1000*10*L1 - P[t]) == (1000/(2*0.0056)*z[t]^2*
L1*(0.316*(1000*z[t]*0.0056/0.00089))^(-0.25)) + 500*01.7*z[t]^2,
(y[t] - z[t])*0.0001 == v'[t], P[t] == p0*v0/v[t]}
/. t -> 0 /. {P -> 100000, v -> 0.001}

and got the following output:

(* {95700. == 5749.49 y^1.75 + 850. y^2,
-50000. == 66854.5 z^1.75 + 850. z^2,
0.0001 (y - z) == v', True} *)

Now let's see if we can find y and z that satisfy the first two algebraic equations:

Solve[95700. == 5749.487805231354 y^1.75 + 850. y^2, y]
(* {{y -> 4.46244}} -- good *)

Solve[-50000. == 66854.50936315529 z^1.75 + 850. z^2, z]
(* {{z -> -0.184864 - 0.821224 I}, {z -> -0.184864 + 0.821224 I}}
-- two complex solutions = bad? *)

Without knowing anything about the physical nature of your equations, I tried flipping the sign of the left hand side of the second equation. That resolved the initial conditions and NDSolve ran. So I'd check to make sure the numbers are right and that the initial conditions can be met.

There are 64 branches of the solution of the system of equations, not all of them are real, four of them are shown in this example

L1 = 5.0;
L2 = 0.43;
d = 0.0056;
g = 9.8;
h1 = 5.0;
h2 = 0.43;
A1 = 0.00001;
A2 = 0.001;
p0 = 100000;
v0 = 0.001;
P[t_] := p0*v0/v[t]

y1[t_, i_] :=
y[t] /. Solve[
P[t] - 1*1000*10*
L2 == (1000/(2*0.0056)*y[t]^2*
L2*(0.316*(1000*y[t]*0.0056/0.00089))^(-0.25)) +
500*01.7*y[t]^2, y[t]][[i]]
z1[t_, i_] :=
z[t] /. Solve[(1000*10*L1 -
P[t]) == (1000/(2*0.0056)*z[t]^2*
L1*(0.316*(1000*z[t]*0.0056/0.00089))^(-0.25)) +
500*01.7*z[t]^2, z[t]][[i]]
sol11 = NDSolveValue[{(y1[t, 1] - z1[t, 1])*0.0001 == v'[t],
v == 0.001}, v[t], {t, 0, 1}];
sol12 = NDSolveValue[{(y1[t, 1] - z1[t, 2])*0.0001 == v'[t],
v == 0.001}, v[t], {t, 0, 1}];
sol13 = NDSolveValue[{(y1[t, 1] - z1[t, 3])*0.0001 == v'[t],
v == 0.001}, v[t], {t, 0, 1}];
sol88 = NDSolveValue[{(y1[t, 8] - z1[t, 8])*0.0001 == v'[t],
v == 0.001}, v[t], {t, 0, 1}];
{Plot[sol11, {t, 0, 1}], Plot[sol12, {t, 0, 1}],
Plot[Re[sol13], {t, 0, 1}], Plot[sol88, {t, 0, 1}]} 