0
$\begingroup$

Here is an ODE solved by ParametricNDSolveValue solver.I got the solution. Now, I need to solve these two equations:

p1 = 10 h - l;
p2 = 10 h1; 

where h and h1 go from 2 to 5. Left to me. l is obtained from solving s1, where

l = 4 π j NIntegrate[r^2 Exp[-First[s1[sol]][r]], {r, 0, R}];

So, I need to get p1 for different values of l. To do that I need to vary z an g; let's say {z, 100, 400, 50} and {g, 0.023, 0.134}. Then for these values, there is a corresponding l to pass to p1. After getting p1 and p2, I want to know where p1 - p2 == 0 by means of inspecting a 3D contour plot. See the following link 1.

Z = 800; 
g = Rationalize[0.0238, 0];
k2 = Rationalize[0.000194519, 0];
ϵ = 10^-4; 
R = Rationalize[1.5472, 0];
j = .00001271;

s1 = 
  ParametricNDSolveValue[
    {y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], 
     y[ϵ] == y0, y'[ϵ] == 0, 
     WhenEvent[r == 1, y'[r] -> y'[r] + Z g]}, 
    {y, y'}, {r, ϵ, R}, {y0}, 
    Method -> "StiffnessSwitching", WorkingPrecision -> 30];

sol = FindRoot[Last[s1[y0]][R], {y0, -11.25, -11.}, 
   Evaluated -> False][[1, 2]]

Plot[First[s1[sol]][r], {r, ϵ, R}, 
  AxesLabel -> {r, y}, 
  ImageSize -> Large, 
  LabelStyle -> {Black, Bold, Medium}]

Plot[
  Evaluate @ Table[TooltipFirst[s1[sol]][r], y0] 
  {z, 100, 400, 50}, {g, 0.023, 0.134}, {r, 0, R}] 

l = 4 π j NIntegrate[r^2 Exp[-First[s1[sol]][r]], {r, 0, R}];
p1 = 10 h - l1 ;
p2 = 10 h1 ; 

ContourPlot3D[
  {p1 == 0, p2 == 0}, {h, -2, 2}, {h1, -2, 2}, {l, -2, 2},  
  MeshFunctions -> {Function[{h, h1, l, f}, p1 - p2]},  
  MeshStyle -> {{Thick, Blue}}, 
  Mesh -> {{0}},  
  ContourStyle ->  Directive[Orange, Opacity[0.5], Specularity[White, 30]]]

My issues that l depends on the s1 and therefore on z, g, etc. I believe I need to use Do or Table, but I do not know how? Please help.

$\endgroup$
  • $\begingroup$ sol has 2 parameters: u0 and c. Therefore l also should depend on u0 and c. And you try to ContourPlot with two additional variables h and h1, I dont't understand what is your goal. $\endgroup$ – Alx Oct 31 '19 at 14:09
  • $\begingroup$ @Alx I updated the Q. $\endgroup$ – all Oct 31 '19 at 17:53
  • $\begingroup$ You write log[2]. Do you really mean Log[2]? $\endgroup$ – m_goldberg Nov 1 '19 at 3:59
  • $\begingroup$ yes I meant 'Log[2]' $\endgroup$ – all Nov 1 '19 at 4:14
  • $\begingroup$ @m_goldberg To make the code simple : assume p1 = 10h-l and p2 = 10 h1. $\endgroup$ – all Nov 1 '19 at 4:38
1
$\begingroup$

I'm not sure I understood your task, but my thought on this.

z = 800; g = Rationalize[0.0238, 0];
k2 = Rationalize[0.000194519, 0];
ϵ = 10^-4; R = Rationalize[1.5472, 0];
j = .00001271;

s1 = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], 
    y[ϵ] == y0, y'[ϵ] == 0, 
    WhenEvent[r == 1, y'[r] -> y'[r] + z1 g1]}, {y, 
    y'}, {r, ϵ, R}, {y0, z1, g1}, 
   Method -> "StiffnessSwitching", WorkingPrecision -> 30];

So, I added two more parameters z1 and g1 (renamed to not to confuse with some initial constants z and g).

Then I define helper function for use inside FindRoot, it accepts numerical arguments only:

s2[y0_?NumericQ, z1_?NumericQ, g1_?NumericQ] := Last[s1[y0, z1, g1]]

Now we can compute y0 values for different z1 and g1:

tab = Flatten[
  Quiet@Table[{y0 /. FindRoot[s2[y0, z, g][R], {y0, -10}], z, g},
  {z, 100, 400, 50}, {g, {0.023, 0.134}}], 1]

and corresponding l values:

l = 4 π j NIntegrate[r^2 Exp[-First[s1[##]][r]], {r, 0, R}] & @@@ tab

What to do next, how to visualize? One can make ContourPlot:

ContourPlot[
 Evaluate@Table[10 h - l1 == 10 h1, {l1, l}], {h, 2, 5}, {h1, 2, 5},
 FrameLabel -> {h, h1}]

enter image description here

Or ListPlot3D:

ListPlot3D[
 Flatten[Table[{h1, l1, h1 + l1/10}, {h1, 2, 5, 0.1}, {l1, l}], 1], 
 AxesLabel -> {h1, l1, h}]

enter image description here

EDIT

To answer additional request: z = 160/h, g = 0.0714/h1. Make corresponding changes:

s1 = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], 
        y[ϵ] == y0, y'[ϵ] == 0, 
        WhenEvent[r == 1, y'[r] -> y'[r] + 160/h 0.0714/h1]}, {y, 
        y'}, {r, ϵ, R}, {y0, h, h1}, 
       Method -> "StiffnessSwitching", WorkingPrecision -> 30];

s2[y0_?NumericQ, h_?NumericQ, h1_?NumericQ] := Last[s1[y0, h, h1]]
ff[h_?NumericQ, h1_?NumericQ] := y0 /. FindRoot[s2[y0, h, h1][R], {y0, -10}]
ll[h_?NumericQ, h1_?NumericQ] := 4 π j NIntegrate[r^2 Exp[-First[s1[ff[h, h1], h, h1]][r]],{r, 0, R}]

Now we are ready to plot contour of 10 h - ll[h, h1] - 10 h1 == 0:

ListContourPlot[
 Flatten[Table[{h, h1, 10 h - ll[h, h1] - 10 h1}, {h, 2, 5, 0.1}, {h1,
     2, 5, 0.1}], 1], Contours -> {0}, ContourShading -> None, 
 ContourStyle -> Thick, FrameLabel -> {h, h1}, BaseStyle -> 12]

enter image description here

FINAL EDIT

k2 = Rationalize[0.000194519, 0];
ϵ = 10^-4; R = Rationalize[1.5472, 0];
j = .00001271;

s1 = ParametricNDSolveValue[{y''[r] + 2 y'[r]/r == k2 Sinh[y[r]], 
    y[ϵ] == y0, y'[ϵ] == 0, 
    WhenEvent[r == 1, y'[r] -> y'[r] + 160/h 0.0714/h1]}, {y, 
    y'}, {r, ϵ, R}, {y0, h, h1}, 
   Method -> "StiffnessSwitching", WorkingPrecision -> 30];

s2[y0_?NumericQ, h_?NumericQ, h1_?NumericQ] := Last[s1[y0, h, h1]]
ff[h_?NumericQ, h1_?NumericQ] := 
 y0 /. FindRoot[s2[y0, h, h1][R], {y0, -10}]
ll[h_?NumericQ, h1_?NumericQ] := 
 4 π j NIntegrate[
   r^2 Exp[-First[s1[ff[h, h1], h, h1]][r]], {r, 0, R}]

To speed-up computation (FindRoot and NIntegrate are very expensive for continuous plotting) I suggest using Interpolation of ll:

llinterp = 
 Interpolation[
  Flatten[Table[{{h, h1}, ll[h, h1]}, {h, 2, 5, 0.1}, {h1, 2, 5, 
     0.1}], 1]]

Now we can define p1 and p2 via llinterp (as of last comment of OP):

p1 = 10 h - llinterp[h, h1] + u
p2 = -u + 10 h1 llinterp[h, h1]

and see what are ContourPlot3Ds:

ContourPlot3D[p1, {h, 2, 5}, {h1, 2, 5}, {u, 2, 10}, Contours -> 5, 
PlotLegends -> Automatic, AxesLabel -> {h, h1, u}]

enter image description here

ContourPlot3D[p2, {h, 2, 5}, {h1, 2, 5}, {u, 2, 10}, Contours -> 5, 
PlotLegends -> Automatic, AxesLabel -> {h, h1, u}]

enter image description here

Colors in lenends correspond to values of p1, p2.

Apparently, with given ranges of variables (h, h1, u) we can not fulfil conditions p1 = 0 and p2 = 0. Moreover, these conditions are conflicting: to decreese p1 one has to decrease h, but for p2 situation is opposite. So, OP should think about conditions and/or variables' ranges.

We can also plot intersection (p1 == p2):

ContourPlot3D[p1 == p2, {h, 2, 5}, {h1, 2, 5}, {u, 2, 10}, 
AxesLabel -> {h, h1, u}]

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks @Alx.I do not see when p1-p2=0 $\endgroup$ – all Nov 1 '19 at 5:45
  • $\begingroup$ Lines in ContourPlot correspond to p1-p2 == 0 which is the same as 10 h - l1 == 10 h1. In ListPlot3D I plot h vs h1 and l1 as a solution to 10 h - l1 == 10 h1 ->h==h1 - l1/10. $\endgroup$ – Alx Nov 1 '19 at 7:04
  • $\begingroup$ I was expecting to see intersection of two surfaces !! which indicates where the zero pint of p1- p2. $\endgroup$ – all Nov 1 '19 at 10:17
  • $\begingroup$ Sorry, don't understand. Let's clarify. p1 and p2 are actually planes in terms of h and h1. You can make this plot: ContourPlot3D[Evaluate[(10 h == 10 h1 - #) & /@ l], {h, 2, 5}, {h1, 2, 5}, {z, -10, 10}] and see that p1 - p2 == 0 generates series of planes, and in cross-section (in h - h1 plane) this gives exactly what I showed in ContourPlot in answer. If it is not what you assumed to be, sorry I have no more ideas. $\endgroup$ – Alx Nov 1 '19 at 12:49
  • $\begingroup$ Thanks Alx. I got your points. What if z and g function of h and h1 : z=160/h and g=0.0714/h1`. Could you please help me with that? $\endgroup$ – all Nov 2 '19 at 0:52

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.