4
$\begingroup$
sol = NDSolve[{y''[x] + (y[x])^2 == 1, y[0] == 2, y'[0] == 3}, y[x], {x, 0, 30}]
g = Plot[Evaluate[{y[x], y'[x], y''[x]} /. sol], {x, 0, 30}, PlotStyle -> Automatic]

I'm trying to coding this equation $y''+ y^2 =1$ with initial value $y(0)= 2$, and $y'(0)= 3$ using NDSolve. Anyone could tell me what is the problem above in the code? Thank you.

$\endgroup$
  • $\begingroup$ It does solve the system, but I get a message which says that there is a singularity at x = 3.95. $\endgroup$ – Lotus Oct 12 '18 at 5:13
  • $\begingroup$ Also, by changing the initial condition to y'[0] == 0 gives a nice periodic solution. Is this what you want ? $\endgroup$ – Lotus Oct 12 '18 at 5:17
  • $\begingroup$ I need it to give me a graph, but I tried old method I learned before, and it didn't work. $\endgroup$ – Mousa Oct 12 '18 at 5:22
  • $\begingroup$ How to change the initial value for y ' [0]=3 to y ' [0]=0 when the given initial in the problem is this one "y ' [0]=3" $\endgroup$ – Mousa Oct 12 '18 at 5:25
  • $\begingroup$ I am not saying that we should change it. I am just pointing out that with y'[0] = 3 we get stiffness issues, which goes away with y'[0] = 0. Understanding of the physics of the problem may explain this behaviour. $\endgroup$ – Lotus Oct 12 '18 at 5:33
5
$\begingroup$

The desired solution can be obtained semi-analytically.

s = DSolveValue[y''[x] + (y[x])^2 == 1, y[x], x] /. {C[1] -> c1, C[2] -> c2}
(* -(-1)^(2/3) 6^(1/3) 
   WeierstrassP[(-(1/6))^(1/3) (c1 + x), {2 (-1)^(2/3) 6^(1/3), c2}] *)

Next, determine the constants by applying the boundary conditions.

bc1 = s /. x -> 0;
bc2 = D[s, x] /. x -> 0;
sc = FindRoot[{bc1 == 2, bc2 == 3}, {c1, I}, {c2, 1}] // Chop
(* {c1 -> -0.523953 + 1.79242 I, c2 -> -10.3333} *)

and plot.

Plot[Chop[s /. sc], {x, 0, 10}, ImageSize -> Large]

enter image description here

The solution is periodic in x. In light of the singularities in the solution, it is not surprising that NDSolve failed for x near 4.

Addendum

The period can be increased by decreasing y'[0] and becomes infinite for y'[0] == 0. A WorkingPrecision of at least 45 is needed in that case.

$\endgroup$
2
$\begingroup$

May not be exactly what OP wants but still:

sol = NDSolve[{y''[x] + (y[x])^2 == 1, y[0] == 2, y'[0] == 0}, 
  y, {x, 0, 30}]

g = Plot[Evaluate[{y[x], y'[x], y''[x]} /. sol], {x, 0, 30}, 
  PlotStyle -> Automatic]

Gives the graphs as required. With y'[0] == 3 as given by OP there seems to be a singularity/stiffness at x = 3.95

$\endgroup$
2
$\begingroup$

just a supplement to @Lotus answer:

Using ParametricNDSolveValue

Y = ParametricNDSolveValue[{y''[x] + (y[x])^2 == 1, y[0] == 2,y'[0] == ys}, y, {x, 0, 30}, ys]

you could do some parameter studies varying ys=y'[0]

Show[{
Plot[Table[Y[p][x], {p, .1, 3. , .1}] , {x, 0, 25},PlotRange -> {-5, 5}] ,
Plot[ Y[0][x] , {x, 0, 25}, PlotStyle -> {Thickness[.01], Red}], 
Plot[ Y[3][x] , {x, 0, 25}, PlotStyle -> { Green },PlotRange -> {-5, 5}]
}]

enter image description here

which show, that only y'[0]==0 gives convergent solution! (The green solution is y'[0]==3 OP asked for )

$\endgroup$
  • $\begingroup$ Increasing WorkingPrecision will give a different plot for y'[0]==0. $\endgroup$ – bbgodfrey Oct 12 '18 at 15:34
  • $\begingroup$ @bbgodfrey Thanks, trying ParametricNDSolveValue[... , WorkingPrecision -> 50] Plot[Y[3][x], {x, 0, 10}, PlotStyle -> {Green}, PlotRange -> {-5, 5} ] doesn't finish ... $\endgroup$ – Ulrich Neumann Oct 12 '18 at 15:54
  • $\begingroup$ Works fine for me with NDSolveValue and Mathematica 11.3. Solution is flat at about -1 for large x, which agrees with the symbolic solution. $\endgroup$ – bbgodfrey Oct 12 '18 at 16:05
  • $\begingroup$ @bbgodfrey Thanks, MMA 11.0.1. NDSolve gives the solution in the range 0<x<3.95 $\endgroup$ – Ulrich Neumann Oct 12 '18 at 19:25

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.