# Using NDSolve for three variables as function of time

I want to be able to plot z, and r. Should I be using DSolve, NDSolve or something else? I've tried to clear it up a bit

Thanks,

Will

ClearAll["Global*"];

(*Constants*)
rho = 1000;
g = 9.81;
Patm = 100000;
gma = 0.0728;
mu = 0.0001;
rs = 287;
T = 293.15;
Psat = 2399;

(*Variables*)
r0 = 25*(10^-6)

(*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 , r == 0.000025 && z == 1 }

sol = NDSolve[eqn, {r[t], Pb[t], z[t]}, t, {Z[t], 1, 0}]

• I am really confused looking at the above code. You have only one differential equation there. dz/dt so why are you trying to solve for 3? You are also defining Pb[t] in terms of r[t] and also defining r[t] as function of Pb[t] which looks like a circle to me. – Nasser Feb 16 '14 at 12:36
• I would like to solve for that one differential equation. But within that I need to solve again for Pb and r? How would i do that? – William Moss Feb 16 '14 at 13:03
• Pl in eqn3 is undefined – qwerty Feb 16 '14 at 13:06
• Yes I have defined Pl now. Thanks. Would it help if i describe the physical system? – William Moss Feb 16 '14 at 13:06
• First, I think this is the call you intend: sol = NDSolve[eqn, {r[t], Pb[t], z[t]}, {t, 1, 0}]. If not, please explain. Second, if that is the call, then Solve[#, Pb] & /@ ({eqn1, eqn2} /. {t -> 0} /. {z -> 1, r -> 0.000025}) shows that your initial conditions are inconsistent. – Michael E2 Feb 16 '14 at 13:52

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.000025, in particular, I suppose for Pb. The solution is to specify both r and Pb.

## 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 == 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 and Pb. 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 -> 1} /. paramSettings,
{r, Pb}, 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.0000247373, Pb == 115696.}
*)


Using both equations in icsOP helps NDSolve initialize the algebraic constraints. We can also see that r == 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.

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 even if we give the solution for r it found above, although it will if we specify the domain Reals.

NSolve[{eqn1, eqn2} /. t -> 0 /.
r -> 0.000024737308346837383 /. {z -> 1} /. paramSettings,
Pb]
(*
{}
*)

• Thanks so much for the help. This is to model bubbles in a column of water. I now just need to account for diffusivity. THANKS AGAIN! – William Moss Feb 17 '14 at 19:35