Skip to main content

plotting - Help in finding a good guesses for ode


This problem solved in 1 by @bbgodfrey. However, when I changed the parameters, the ode becomes very stiff and I did not know which best guess I need to use for y0. I am looking for :


 ParallelTable[Quiet[p[Z, g, k2, 1.5]], {g, 0.01, 0.03, .02}, {Z, 1000, 4000, 500},{k2, 0.0002, 100, 50}]

Does anyone know a method to estimate a good initial guess? I do not want to just put a random tries every time. I tried -.1 < y0 < -.08,but did not find a solution.



Here is my code:


p[Z0_, g0_, k0_, R0_] := Block[{Z = Z0, 
      g = Rationalize[g0, 0], k2 = Rationalize[k0, 0], ϵ = 10^-4, R = Rationalize[R0, 0]}, 
      ps = 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, ϵ], R}, {y0}, 
        Method -> "StiffnessSwitching", WorkingPrecision -> 20];
      sol = FindRoot[Last[ps[y0]], {y0, -1}, Evaluated -> False][[1, 2]];
      tot = 4 π sol NIntegrate[r^2 Exp[-First[ps[sol]][r]], {r, ϵ, R}];
      L = Z/g*tot]
DistributeDefinitions[p];

tab = ParallelTable[Quiet[p[Z, g, 0.0002, 1.5]], {Z, 500, 800, 50}, {g, 0.02, 0.06, .005}];
ListContourPlot[tab, DataRange -> {{.02, .06}, {500, 800}}, PlotLegends -> Automatic, 
    FrameLabel -> {g, Z}, ImageSize -> Large, LabelStyle -> {Bold, Black, 15}]

Answer



Simply stated, a "good" initial guess for y0 in FindRoot is one for which ParametricNDSolve can integrate the ODE in the question from r = ϵ to r = R without errors. For the parameter ranges used in question 210296, -1 or any other negative number roughly in the range {-5, 0} worked well. The present question, however, introduces a much wider set of parameters, some of which accommodate only tiny ranges for the initial guesses. However, if the desired guess can be bracketed in y0, then the bracketed region can systematically be shrunk until a valid guess is obtained.


Specifically, fy in the code block below begins by bifurcating, yc = (yl + yu)/2, the bracketing region, and checks by means of tst = First[ps[yc]]["Domain"][[1, 2]] to determine whether ParametricNDSolve is able to integrate the ODE all the way to R for y0 = yc. If so, fy has found a good y00, which it Throws to p and terminates. Otherwise, it takes advantage of the fact that y'[tst] is negative when yc is more negative than the band of good guesses and is positive when yc is more positive than the band of good guesses in order to reduce the size of the bracketing region by a factor of two, after which the process is iterated until a good guess is found. Typically, less than ten iterations are required to find y00.


p[Z0_, g0_, k0_, R0_] := Block[{Z = Z0, g = Rationalize[g0, 0], k2 = Rationalize[k0, 0], 
        ϵ = 10^-4, R = Rationalize[R0, 0], yl = -8, yu = 0, ps, fy, y00, sol}, 
    ps = 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 -> 20];
    fy := (Do[yc = (yl + yu)/2; 
        tst = First[ps[yc]]["Domain"][[1, 2]];
        If[tst == R, Throw[yc]]; 
        If[Last[ps[yc]][tst] > 0, yu = yc, yl = yc], {i, 50}]; yc);
    y00 = Which[R == First[ps[yl]]["Domain"][[1, 2]], yl, 
        R == First[ps[yu]]["Domain"][[1, 2]], yu, True, Catch[fy]];
    sol = FindRoot[Last[ps[y0]][R], {y0, y00}, Evaluated -> False][[1, 2]];
    L = Z/g 4 π sol NIntegrate[r^2 Exp[-First[ps[sol]][r]], {r, ϵ, R}]]


Applying p to the parameter set introduced by the present question yields.


DistributeDefinitions[p];
ParallelTable[Quiet[p[Z, g, k2, 1.5]], {Z, 1000, 4000, 500}, {g, 0.01, 0.03, .02}, 
    {k2, 0.0002, 100, 50}]

(* {{{-1.36776*10^12, -14119.8}, {-1.422*10^12, -16956.}}, 
    {{-3.13388*10^12, -33652.5}, {-3.2242*10^12, -40906.5}}, 
    {{-5.62603*10^12, -63057.8}, {-5.79668*10^12, -75910.8}}, 
    {{-8.84774*10^12, -102510.}, {-9.06488*10^12, -122341.}}, 
    {{-1.2798*10^13, -152604.}, {-9.58419*10^10, -179866.}}, 

    {{-1.74782*10^13, -213431.}, {-6.31954*10^10, -248535.}}, 
    {{-2.2886*10^13, -287387.}, {-5.70976*10^10, -327845.}}} *)

The computation requires about ninety seconds on my six-processor computer, with only a small portion of that time used to obtain initial guesses. Note that the iteration over Z, which has most elements, is placed at the leftmost position to maximize the advantage of parallelization. Placing it as in the question increases the runtime to four minutes.


In closing, I would make two points. First, the "guess finder function", fy, fails if {yl, yu} does not bracket at least part of the band of good guesses. If fy does not converge to a good guess, sufficiently decreasing yl should solve the problem. Second, WorkingPrecision may need to be increased for some parameters to improve accuracy.


Comments

Post a Comment

Popular posts from this blog

plotting - How to draw lines between specified dots on ListPlot?

I would like to create a plot where I have unconnected dots and some connected. So far, I have figured out how to draw the dots. My code is the following: ListPlot[{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4,13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full] I have thought using ListLinePlot command, but I don't know how to specify to the command to draw only selected lines between the dots. Do have any suggestions/hints on how to do that? Thank you. Answer One possibility would be to use Epilog with Line : ListPlot[ {{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4, 13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full, Epilog -> { Line[ ...
Post a Comment
Read more

equation solving - Invert and fit implicitly defined curve

I need to fit an implicitly defined curve. I thought I could get some data out of Solve , and then using FindFit . Therefore, I would like to find the relation the parametric curve defined by $F(x,y)=0$: Solve[-(1/2) + 1/2 (0.41202 BesselK[0, 0.1 Sqrt[x^2 + y^2]] + (0.101483 x BesselK[1, 0.1 Sqrt[x^2 + y^2]])/Sqrt[x^2 + y^2]) == 0, y] But I can't get an output: Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help. >> Edit: In particular, I would like to fit the data coming from the curve with the expression of another curve, and not with a function $f(x)$. In particular, since this clearly looks like a cardioid , I would like it to fit to something like it. What other strategies could I try?
Post a Comment
Read more

dynamic - How can I make a clickable ArrayPlot that returns input?

I would like to create a dynamic ArrayPlot so that the rectangles, when clicked, provide the input. Can I use ArrayPlot for this? Or is there something else I should have to use? Answer ArrayPlot is much more than just a simple array like Grid : it represents a ranged 2D dataset, and its visualization can be finetuned by options like DataReversed and DataRange . These features make it quite complicated to reproduce the same layout and order with Grid . Here I offer AnnotatedArrayPlot which comes in handy when your dataset is more than just a flat 2D array. The dynamic interface allows highlighting individual cells and possibly interacting with them. AnnotatedArrayPlot works the same way as ArrayPlot and accepts the same options plus Enabled , HighlightCoordinates , HighlightStyle and HighlightElementFunction . data = {{Missing["HasSomeMoreData"], GrayLevel[ 1], {RGBColor[0, 1, 1], RGBColor[0, 0, 1], GrayLevel[1]}, RGBColor[0, 1, 0]}, {GrayLevel[0], GrayLevel...
Post a Comment
Read more