Once the trivial things are fixed -- I borrowed the definition of Pl
from the Wolfram Community posting -- the problem is that NDSolve
has trouble setting up the initial conditions from r[0] == 0.000025
, in particular, I suppose for Pb[0]
. The solution is to specify both r[0]
and Pb[0]
.
OP's setup
I converted each parameter setting to a Rule
, so I could play with them more easily.
(*equations*)
Pl = rho*g*z[t] + Patm;
eqn1 = Pb[t] == rho*g*z[t] + Patm + (2*gma)/r[t];
eqn2 = r[t] == (((1 + (Psat/Pb[t]))*(rho*g*z[t] + Patm)*r0^3)/Pb[t])^(1/3);
eqn3 = z'[t] == -(2/9)*(1/mu)*(r[t]^2)*(Pl - Pb[t]/(rs*T));
eqn = {eqn1 && eqn2 && eqn3, z[0] == 1};
params12 = Union@Cases[{eqn1, eqn2}, _Symbol, Infinity] /. t -> Sequence[];
params = Union@Cases[{eqn1, eqn2, eqn3}, _Symbol, Infinity] /. t -> Sequence[];
paramSettings = Hold[
rho = 1000;
g = 9.81;
Patm = 100000;
gma = 0.0728;
mu = 0.0001;
rs = 287;
T = 293.15;
Psat = 2399;
r0 = 25*(10^-6)] /. Set -> Rule /. CompoundExpression -> List //
ReleaseHold
(*
{rho -> 1000, g -> 9.81, Patm -> 100000, gma -> 0.0728, mu -> 0.0001,
rs -> 287, T -> 293.15, Psat -> 2399, r0 -> 1/40000}
*)
Solving the OP's setup
We solve the algebraic constraints at t -> 0
for the initial values of r[0]
and Pb[0]
. The warning makes me worry about the sensitivity of the results to the uncertainty of the parameter settings. The warning does not suggest that, but converting to exact coefficients could hide it.
icsOP = Equal @@@
First@NSolve[{eqn1, eqn2} /. t -> 0 /. {z[0] -> 1} /. paramSettings,
{r[0], Pb[0]}, Reals]
NSolve::ratnz: NSolve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>
(*
{r[0] == 0.0000247373, Pb[0] == 115696.}
*)
Using both equations in icsOP
helps NDSolve
initialize the algebraic constraints. We can also see that r[0] == 0.000025
is a little off, which is why NDSolve
won't solve the system with the original setup. It still seems to have trouble if I give it just the one for r[0]
.
solOP = First @ NDSolve[{eqn, icsOP} /. paramSettings, {r[t], Pb[t], z[t]}, {t, 1, 0}]
(*
{r[t] -> InterpolatingFunction[{{0., 1.}}, <>][t],
Pb[t] -> InterpolatingFunction[{{0., 1.}}, <>][t],
z[t] -> InterpolatingFunction[{{0., 1.}}, <>][t]}
*)
Below we see that NSolve
does not find a solution for Pb[0]
even if we give the solution for r[0]
it found above, although it will if we specify the domain Reals
.
NSolve[{eqn1, eqn2} /. t -> 0 /.
r[0] -> 0.000024737308346837383` /. {z[0] -> 1} /. paramSettings,
Pb[0]]
(*
{}
*)
dz/dt
so why are you trying to solve for 3? You are also definingPb[t]
in terms ofr[t]
and also definingr[t]
as function ofPb[t]
which looks like a circle to me. $\endgroup$sol = NDSolve[eqn, {r[t], Pb[t], z[t]}, {t, 1, 0}]
. If not, please explain. Second, if that is the call, thenSolve[#, Pb[0]] & /@ ({eqn1, eqn2} /. {t -> 0} /. {z[0] -> 1, r[0] -> 0.000025})
shows that your initial conditions are inconsistent. $\endgroup$