0
$\begingroup$

Below is an ODE with BC define as x[R]=0, x'[0]=0, x'[R]=0 and parameter n. The ODE is stiff at certain data, and I need to see the behavior of x' and x for given parameters. I appreciate any help.

c = -1;
r1 = 0.8;
r2 = 1;
R = 1.29;
f[r_] := Piecewise[{{0, 0 <= r <= r1}, {900/(1 - r1^3), 
    r1 < r <= 1}, {0, 1 < r <= R}}]


ps = ParametricNDSolveValue[{x''[r] + (1/r) x'[r] == 
    c n Exp[-x[r]] + f[r], x'[0] == 0, x[0] == x0}, {x, x'}, {r, 0, 
   R}, {x0,n}, Method -> "StiffnessSwitching"]

ff = FindRoot[{Last[ps[x0,n]][R] == 0, First[ps[x0,n]][R] == 0}, {x0, -2}]
$\endgroup$
  • $\begingroup$ I have not yet found a convergent solution. It is possible that it is necessary to change the boundary conditions. $\endgroup$ – Alex Trounev Oct 10 '18 at 18:10
  • $\begingroup$ @AlexTrounev see my update $\endgroup$ – user60416 Oct 10 '18 at 19:09
  • $\begingroup$ The figure shows R>r2, and in the data otherwise r2=125>R=1.29. $\endgroup$ – Alex Trounev Oct 11 '18 at 3:53
  • $\begingroup$ Here r1/r2=0.8 and R/r2 =1.29, so instead the integral from r1<r<r2 is it from r1<r<1 where 1 = r2. $\endgroup$ – user60416 Oct 11 '18 at 4:40
  • $\begingroup$ What is the meaning of the parameter с, if the solution of the equation depends on the product c * n? $\endgroup$ – Alex Trounev Oct 11 '18 at 4:47
3
$\begingroup$

The solution of the equation depends on the product c * n, but not on c. Therefore, we can set c = -1, since the parameter n is found in the solution process.

c = -1;
r1 = 8/10;
r2 = 1;
R = 129/100;
r0 = 10^-6;
f[r_] := Piecewise[{{0, 0 <= r <= r1}, {900/(1 - r1^3), 
    r1 < r <= 1}, {0, 1 < r <= R}}]


ps = ParametricNDSolveValue[{x''[r] + (1/r) x'[r] == 
     c n Exp[-x[r]] + f[r], x'[r0] == 0, x[r0] == x0}, 
   x, {r, r0, R}, {x0, n}, Method -> "StiffnessSwitching", 
   WorkingPrecision -> 30];
ps1 = ParametricNDSolveValue[{x''[r] + (1/r) x'[r] == 
    c n Exp[-x[r]] + f[r], x'[r0] == 0, x[r0] == x0}, 
  x', {r, r0, R}, {x0, n}, Method -> "StiffnessSwitching", 
  WorkingPrecision -> 30]

ff = FindRoot[{ps[x0, n][R] == 0, ps1[x0, n][R] == 0}, {x0, 
   1.277171`30}, {n, 393/10}]



{x0 -> 1.27686, n -> 39.2879}

{Plot[Evaluate[ps[Last[ff[[1]]], Last[ff[[2]]]][r]], {r, r0, R}, 
  PlotRange -> All, AxesLabel -> {"r", "x"}], 
 Plot[Evaluate[ps1[Last[ff[[1]]], Last[ff[[2]]]][r]], {r, r0, R}, 
  PlotRange -> All, AxesLabel -> {"r", "x'"}]}

fig1

$\endgroup$
  • $\begingroup$ How do you estimate theses numbers {x0, 1.277171`30}, {n, 393/10}]? $\endgroup$ – user60416 Oct 11 '18 at 19:23
  • $\begingroup$ Would you please explain to me this line Evaluate[ps[Last[ff[[1]]], Last[ff[[2]]]][r]]? like what is 2 refer to? $\endgroup$ – user60416 Oct 12 '18 at 1:13
  • $\begingroup$ To find the initial values {x0,n}, I used an iteration prediction procedure. The second question is not clear to me. Why don't you just see what corresponds to what? ps->x,ps1->x',Last[ff[[1]]]-> 1.27686,Last[ff[[2]]]->39.2879 $\endgroup$ – Alex Trounev Oct 12 '18 at 2:57
  • $\begingroup$ Thanks! Could you share your code of (iteration prediction procedure)? $\endgroup$ – user60416 Oct 12 '18 at 3:31
  • $\begingroup$ Unfortunately I can not give the code. But it is not difficult to write in the Wolfram Language. $\endgroup$ – Alex Trounev Oct 12 '18 at 3:35

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.