1
$\begingroup$

I'm trying to solve a 2nd order related ODE system in two region using Piecewise. BC: 1) y'[0]=0 and at y'[1.6]=0 while at 1 y'[1] is continues.y[1.6]=0andy[0]= y0.

I'm using to use NDSolve as shown below, but the code did not solve the system

How could I write code using NDSolve that the solution match at `y'[1] in the two region? The code is:

 \[Epsilon] = $MachineEpsilon ;
f[x_] := Piecewise[{{2.7, 0 <= x <= 1}, {0, 1 < x <= 1.6}}];
s = NDSolve[{y''[x] + 2 y'[x]/x == - Cosh[-y[x]] + f[x], 
   y[1] == y0, y'[\[Epsilon]] == 0}, y, {x, \[Epsilon], 1.6},
   Method -> "StiffnessSwitching", WorkingPrecision -> 40]
Plot[%, {x, 0, 1.6}, AxesLabel -> {x, y}]
$\endgroup$
  • $\begingroup$ You write y[\[Epsilon]] == y0, but there is no y0 in your code. it will be better to make the code self contained so it it easy to copy and paste it in one copy. $\endgroup$ – Nasser Jun 11 at 4:25
  • $\begingroup$ y0 is the value of function at \[Epsilon]that need to evaluate. $\endgroup$ – user66021 Jun 11 at 4:34
  • $\begingroup$ y0 is the value of function at \[Epsilon]that need to evaluate I am not following. You are giving an initial condition for the dependent variable at some location as unknown? This is NDSolve not DSolve. So I think it should be known value there. $\endgroup$ – Nasser Jun 11 at 4:37
  • $\begingroup$ Ok. ODE has to region, and my goals to solve ode in each region and match the solution at Y(1). So, how could I do this? $\endgroup$ – user66021 Jun 11 at 4:46
  • $\begingroup$ Moreover, I need to pass an initial estimate to the function at x=1, and from that using find root. But I did know how? $\endgroup$ – user66021 Jun 11 at 4:53
2
$\begingroup$
Plot[{Piecewise[{{2.7, 0 <= x <= 1}, {0, 1 < x <= 1.6}}], 
  2.7/(1 + Exp[200 (x - 1)])}, {x, 0, 1.6}]

enter image description here

\[Epsilon] = $MachineEpsilon;
f[x_] := 2.7/(1 + Exp[200 (x - 1)]);
s[c_] := y /. 
   First@NDSolve[{y''[x] + 2 y'[x]/x == c*-Cosh[-y[x]] + f[x], 
      y[1.6] == 0, y'[\[Epsilon]] == 0}, y, {x, \[Epsilon], 1.6}, 
     Method -> "StiffnessSwitching", AccuracyGoal -> 30];
sol = s[1];
Plot[sol[x], {x, 0, 1.6}, AxesLabel -> {x, y}]

enter image description here


With variable parameter c,

Do[
  sols[i] = s[i];, {i, 0, 1, 0.1}];

Plot[MapThread[sols[#][x] &, {Range[0, 1, 0.1]}], {x, 0, 1.6}, 
 AxesLabel -> {x, y}]

Mathematica graphics


With Manipulate

Manipulate[
 Plot[sols[c][x], {x, 0, 1.6}, AxesLabel -> {x, y}], {c, 0, 1, 0.1}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ let say I have a parameter c y''[x] + 2 y'[x]/x == -Cosh[-y[x]] + f[x] > y''[x] + 2 y'[x]/x == -c Cosh[-y[x]] + f[x] ?How to insert it in your NDSolve method ? $\endgroup$ – user66021 Jun 11 at 19:13
  • $\begingroup$ @user66021 would you ask/post it as another(new) question? $\endgroup$ – Xminer Jun 11 at 19:15
  • $\begingroup$ Is it worth to put a new questions since I'm just adding a parameter. ? $\endgroup$ – user66021 Jun 11 at 19:17
  • $\begingroup$ @user66021 ok,I've understand what you mean from this comment. $\endgroup$ – Xminer Jun 11 at 19:41
  • $\begingroup$ Yes, Would you please do it in NDSolve? $\endgroup$ – user66021 Jun 11 at 21:23
2
$\begingroup$

You can solve your problem using ParametricNDSolve.

Because MMA currently (v 11.0.1) isn't able to solve boundary value problems I define an initial value problem at the boundary x=1.6 with unknown slope y'[1.6]==ys1 thereby avoiding the singularity at x==0

Y = ParametricNDSolveValue[{y''[x] + 2 y'[x]/x == -Cosh[-y[x]] + f[x], y[1.6] == 0, y'[1.6] == ys1}, y, {x, 0, 1.6}, {ys1}] 

Now try to find ys1 such that y'[0]->0

opt = NMinimize[{1, Y[ys1]'[0] == 0 , -.5 < ys1 < 0}, {ys1}]
Plot[Y[ys1][x] /. opt[[2]] , {x, 0, 1.6}, Evaluated -> True]

enter image description here

In this approach there is no need to introduce small \[Epsilon], special method (NDSolve), very high WorkingPrecision or simplified f[x]

$\endgroup$
  • $\begingroup$ Thanks for those solutions. What will happen if y''[x] + 2 y'[x]/x == -Cosh[-y[x]] + f[x] has a parameter (c) > y''[x] + 2 y'[x]/x == -c Cosh[-y[x]] + f[x] ?How to insert it in NDSolve method ? $\endgroup$ – user66021 Jun 11 at 8:38
  • $\begingroup$ ParametricNDSolveValue[..., {ys1,c}] $\endgroup$ – Ulrich Neumann Jun 11 at 8:43
  • $\begingroup$ Thanks @Ulrich Neumann. $\endgroup$ – user66021 Jun 11 at 8:46

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.