# 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. Commented Dec 14, 2018 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
Commented Dec 14, 2018 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. Commented Dec 15, 2018 at 0:56
• @MichaelE2 It is Pi, sorry about that. But still I am getting errors .
– JaJ
Commented Dec 15, 2018 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
Commented Dec 15, 2018 at 3:49
• What is this mem never saw this before?
– JaJ
Commented Dec 15, 2018 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
Commented Dec 15, 2018 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. Commented Dec 15, 2018 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}]