13

I am trying to solve the following differential equation-

k = NDSolve[{y'''[x] == -76*Sin[y[x]], y[0] == Pi/4, y'[0] == 0, 
   y[1.85] == 0}, y, {x, 0, 2}, AccuracyGoal -> 10, 
  PrecisionGoal -> 30]

I am getting a solution but when I plot a graph, it does not satisfy the boundary conditions. I am getting the following two errors- FindRoot::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option. >>

NDSolve::berr: There are significant errors {-0.239832,-0.0591905,-0.0015434} in the boundary value residuals. Returning the best solution found. >>

The funny thing is when I reduce the coefficient of Sin[y[x]] it works fine. Can anyone help me with this? Thanks,

Community
  • 1
valmiki
  • 131
  • 3
  • You could start by trying a "shooting method" approach, where you use only initial conditions - i.e. replace y[1.85] == 0 by y''[0] == c, and then vary c to see what value of y[1.85] it gives rise to. I found that y[1.85] rapidly switches between negative and positive values as c increases from 24.4 to 24.5, so your solution lies somewhere in that range. The optimisation of c should be automated, but I haven't got time right now. – Stephen Luttrell Jan 23 '13 at 16:55

1 Answers1

11

You could use the "shooting method", as described in my earlier comment, to find the solution as follows:

Clear[error, y, x, c];

Define a function to minimise the predicted value of y[1.85]^2 if you assume the initial condition y''[0] == c:

error[c_?NumericQ] := 
  y[1.85]^2 /. 
   NDSolve[{y'''[x] == -76*Sin[y[x]], y[0] == Pi/4, y'[0] == 0, 
      y''[0] == c}, y, {x, 0, 2},
      AccuracyGoal -> 10, PrecisionGoal -> 30][[1]];

[EDIT: I really should have used y[1.85] rather than y[1.85]^2, because this is a better match for the FindRoot that I use below, but in this case FindRoot finds the (second order) zero anyway. The squared version y[1.85]^2 is what you would use if you were then going to use NMinimize - i.e. the minimum is not necessarily zero.]

Find the solution given the starting value for c mentioned in my earlier comment - it doesn't have to be this accurate a guess.

c0 = c /. FindRoot[error[c], {c, 24.4}]

(* 24.4354 *)

This tells you that, for the initial conditions y[0] == Pi/4 and y'[0] == 0, the end point condition y[1.85] == 0 is equivalent to the extra initial condition y''[0] == 24.4354.

Plug this solution for c back in

soln = NDSolve[{y'''[x] == -76*Sin[y[x]], y[0] == Pi/4, y'[0] == 0, 
  y''[0] == c0}, y, {x, 0, 2},
  AccuracyGoal -> 10, PrecisionGoal -> 30][[1]];

and plot the result

Plot[y[x] /. soln // Evaluate, {x, 0, 2}]

Mathematica graphics

Stephen Luttrell
  • 5,044
  • 1
  • 19
  • 18