# A complicated use of NDSolve

I am looking for a particular type of plots. OK, let me explain,

Here is my DE (I don't know what to call it), where, $$\frac{\partial p}{\partial r} =0$$, $$r_1=\epsilon$$, $$r_2=1+\phi\cos(2\pi z)$$. We have another quantity, $$Q(z)=Q^*-(1+\phi^2/2)+\epsilon^2$$,

$$Q(z)=2\int_{r1}^{r2}r w[r,z]dr$$,

and finally, we have another quantity,

$$\Delta P=\int_{0}^{1}\frac{\partial p}{\partial z}dz$$. The paper, I am reviewing, does this,

1. First solve the original PDE for $$w[r,z]$$ by treating it as an ODE assuming $$p_z$$ to be a constant.

2. Integrate $$w[r,z]$$ w.r.t $$r$$, then sub the expression for $$Q(z)$$ and then solve for $$p_z$$.

3. Finally, integrate $$p_z$$ w.r.t $$z$$ to get $$\Delta P$$.

The authors find analytical solution, which then make the whole process easy but I am looking for a numerical one.

Here is my trying but I am unable to get anywhere with it.

r2 = 1 + phi*Cos[2*Pi*z]; K1 = 0.2; epsilon = 0.32; phi = 0.4;

sol[P1_, z_] :=
First@NDSolve[{P1 == 1/r*D[r*D[w[r], r], r] - 1/K1*(w[r] + 1),
w[epsilon] == -1, w[1 + phi*Cos[2*Pi*z]] == -1},
w, {r, epsilon, 1 + phi*Cos[2*Pi*z]}]

PP1[z_, Q_] :=
NSolve[(Q - (1 + phi^2/2) + epsilon^2) -
2*NIntegrate[
r*w[r] /. sol[P1, z], {r, epsilon, 1 + phi*Cos[2*Pi*z]}], P1]

Plot[PP1[z, 0.1], {z, 0, 1}]

DP[Q_] := NIntegrate[PP1[z, 0.1], {z, 0, 1}]

Plot[DP[Q], {Q, -1, 1}]

• I think of NSolve as a hybrid symbolic-numeric solver. You probably want FindRoot instead. If I have time later, I'll take a look. – Michael E2 Dec 14 '18 at 15:27
• You should also be careful when calling a minimization inside another minimization for example check this thread : mathematica.stackexchange.com/questions/105948/… – user59583 Dec 14 '18 at 16:02
• Capital Phi is undefined and therefore the numeric functions cannot solve equations with symbolic parameters. Is it meant to be lowercase phi. – Michael E2 Dec 15 '18 at 0:56
• @MichaelE2 It is Pi, sorry about that. But still I am getting errors . – JaJ Dec 15 '18 at 3:01

A thorough use of ?NumericQ as well as a few other adjustments gives a result:

ClearAll[sol, PP1, DP];
r2 = 1 + phi*Cos[2*Pi*z]; K1 = 0.2; epsilon = 0.32; phi = 0.4;

mem : sol[P1_?NumericQ, z_?NumericQ] :=
mem = First@
NDSolve[{P1 == 1/r*D[r*D[w[r], r], r] - 1/K1*(w[r] + 1),
w[epsilon] == -1, w[1 + phi*Cos[2*Pi*z]] == -1},
w, {r, epsilon, 1 + phi*Cos[2*Pi*z]}]

PP1[z_?NumericQ, Q_?NumericQ] := Module[{obj, P1},
obj[P1_?NumericQ] := (Q - (1 + phi^2/2) + epsilon^2) -
2*NIntegrate[
r*w[r] /. sol[P1, z], {r, epsilon, 1 + phi*Cos[2*Pi*z]}];
P1 /. FindRoot[obj[P1], {P1, 1}]
];

DP[Q_?NumericQ] := NIntegrate[PP1[z, Q], {z, 0, 1}]

Plot[PP1[z, 0.1], {z, 0, 1}, MaxRecursion -> 0] Plot[DP[Q], {Q, -1, 1}, MaxRecursion -> 0]
(* takes too long *)

• WoW! thanks. Is it possible to make it faster? I mean plotting DP. – JaJ Dec 15 '18 at 3:49
• What is this mem never saw this before? – JaJ Dec 15 '18 at 3:50
• With Table I can generate data to plot DP. So this problem is solved. Now I am left with another one, how I can plot w[r] vs r? – JaJ Dec 15 '18 at 5:25
• mem : ... is a Pattern name matching the function call and is used for memoization. It stores the result so that it will not be recomputed. It's probably not help with speed at all here, because it's unlikely that the functions are being calculated at the same inputs more than once. I forgot it was there or I would have removed it. -- As for other speed-ups, it's unlikely considering how many times NDSolve is called. – Michael E2 Dec 15 '18 at 13:39

To be sure that the numerical solution coincides with the analytical one, we will find two solutions and compare them. Analytical solution

eq = D[r*D[w[r], r], r]/r - (w[r] + 1)/k - pz == 0;
bc = {w[r1] == -1, w[r2] == -1};
sol = DSolve[{eq, bc}, w[r], r];

Q = 2*Integrate[r*w[r] /. First[sol], {r, r1, r2}];

PZ = pz /.
First[Solve[Q - Qa + 1 + \[CurlyPhi]^2/2 - \[Epsilon]^2 == 0, pz]];

PZ1 = PZ /. {r1 -> \[Epsilon], r2 -> 1 + \[CurlyPhi]*Cos[2*Pi*z]};

pic1 = Table[
Plot[Re[PZ1] /. {\[CurlyPhi] -> .4, \[Epsilon] -> .32, k -> K1,
Qa -> Q1}, {z, 0, 1}, AxesLabel -> {"z", "dp/dz"},
PlotLabel ->
Grid[{{"\!$$\*OverscriptBox[\(Q$$, $$_$$]\)=", Q1}, {"K=",
K1}}]], {Q1, {0.1, .3}}, {K1, {0.05, .25}}]

pic2 = Table[
Plot[Re[PZ1] /. {\[CurlyPhi] -> phi, \[Epsilon] -> eps, k -> .1,
Qa -> .1}, {z, 0, 1}, AxesLabel -> {"z", "dp/dz"},
PlotLabel ->
Grid[{{"\[CurlyPhi]=", phi}, {"\[Epsilon]=",
eps}}]], {phi, {0.3, .4}}, {eps, {0.38, .42}}] Numerical solution has already been obtained by Michael E2, use it and compare two solutions.

ClearAll[sol, PP1, DP];
r2 = 1 + phi*Cos[2*Pi*z]; K1 = 0.2; epsilon = 0.32; phi = 0.4;

sol[P1_?NumericQ, z_?NumericQ] :=
First@NDSolve[{P1 == 1/r*D[r*D[w[r], r], r] - 1/K1*(w[r] + 1),
w[epsilon] == -1, w[1 + phi*Cos[2*Pi*z]] == -1},
w, {r, epsilon, 1 + phi*Cos[2*Pi*z]}]

PP1[z_?NumericQ, Q_?NumericQ] :=
Module[{obj, P1},
obj[P1_?NumericQ] := (Q - (1 + phi^2/2) + epsilon^2) -
2*NIntegrate[
r*w[r] /. sol[P1, z], {r, epsilon, 1 + phi*Cos[2*Pi*z]}];
P1 /. FindRoot[obj[P1], {P1, 1}]];

DP[Q_?NumericQ] := NIntegrate[PP1[z, Q], {z, 0, 1}]

Plot[{PP1[z, 0.1],
Re[PZ1] /. {\[CurlyPhi] -> .4, \[Epsilon] -> .32, k -> .2,
Qa -> .1}}, {z, 0, 1}, MaxRecursion -> 0,
PlotLegends -> {"Numerical", "Analytical"},
PlotStyle -> {Thin, Dashed}] 