13
$\begingroup$

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,

Improve this question
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Stephen Luttrell
    Jan 23, 2013 at 16:55

1 Answer 1

Reset to default
11
$\begingroup$

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

Improve this answer
$\endgroup$
1
  • $\begingroup$ Export["/where/you/want/image.png", %] Or use @Szabolcs automated uploader $\endgroup$
    – ssch
    Jan 23, 2013 at 18:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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