3
$\begingroup$

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),

enter image description here

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}]
Improve this question
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Michael E2 Dec 14 '18 at 15:27
  • $\begingroup$ You should also be careful when calling a minimization inside another minimization for example check this thread : mathematica.stackexchange.com/questions/105948/… $\endgroup$ – user59583 Dec 14 '18 at 16:02
  • $\begingroup$ Capital Phi is undefined and therefore the numeric functions cannot solve equations with symbolic parameters. Is it meant to be lowercase phi. $\endgroup$ – Michael E2 Dec 15 '18 at 0:56
  • $\begingroup$ @MichaelE2 It is Pi, sorry about that. But still I am getting errors . $\endgroup$ – JaJ Dec 15 '18 at 3:01
4
$\begingroup$

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]

enter image description here

Plot[DP[Q], {Q, -1, 1}, MaxRecursion -> 0]
(* takes too long *)
Improve this answer
$\endgroup$
  • $\begingroup$ WoW! thanks. Is it possible to make it faster? I mean plotting DP. $\endgroup$ – JaJ Dec 15 '18 at 3:49
  • $\begingroup$ What is this mem never saw this before? $\endgroup$ – JaJ Dec 15 '18 at 3:50
  • $\begingroup$ 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? $\endgroup$ – JaJ Dec 15 '18 at 5:25
  • $\begingroup$ 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. $\endgroup$ – Michael E2 Dec 15 '18 at 13:39
2
$\begingroup$

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}}]

fig1 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}]

fig2

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.