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}]
  • I have not yet found a convergent solution. It is possible that it is necessary to change the boundary conditions. – Alex Trounev Oct 10 at 18:10
  • @AlexTrounev see my update – user60416 Oct 10 at 19:09
  • The figure shows R>r2, and in the data otherwise r2=125>R=1.29. – Alex Trounev Oct 11 at 3:53
  • 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. – user60416 Oct 11 at 4:40
  • What is the meaning of the parameter с, if the solution of the equation depends on the product c * n? – Alex Trounev Oct 11 at 4:47

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

  • How do you estimate theses numbers {x0, 1.277171`30}, {n, 393/10}]? – user60416 Oct 11 at 19:23
  • Would you please explain to me this line Evaluate[ps[Last[ff[[1]]], Last[ff[[2]]]][r]]? like what is 2 refer to? – user60416 Oct 12 at 1:13
  • 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 – Alex Trounev Oct 12 at 2:57
  • Thanks! Could you share your code of (iteration prediction procedure)? – user60416 Oct 12 at 3:31
  • Unfortunately I can not give the code. But it is not difficult to write in the Wolfram Language. – Alex Trounev Oct 12 at 3:35

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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