Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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,

share|improve this question
    
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
add comment

1 Answer

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

share|improve this answer
    
Export["/where/you/want/image.png", %] Or use @Szabolcs automated uploader –  ssch Jan 23 '13 at 18:03
    
Thanks Stephen, –  valmiki Jan 24 '13 at 7:55
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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