Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I want to solve two coupled equations with NDSolve,

t x'[t] == -x[t] + y[t], t y'[t] == -5 t^2/x[t]^2 + x[t] - y[t], 
x[0] == y[0], x[1] == 1

It's easy to see that at t = 0 the derivative is undetermined, so NDSolve fails.

The Mathematica help center told me one possibility is to start at a small ε > 0 instead of 0, and another possibility is to use the option SolveDelayed -> True to avoid the singularity in the solved form of the equations.

Unfortunately neither of the methods worked for me. What should I do?

share|improve this question
up vote 11 down vote accepted

You can try "shooting method":

eqn1 = t x'[t] - (-x[t] + y[t]);
eqn2 = t y'[t] - (-5 t^2/x[t]^2 + x[t] - y[t]);
sol = NDSolve[{eqn1 == 0, eqn2 == 0, x[0] == y[0], x[1] == 1}, {x, y}, {t, 0, 1}, 
              Method -> {"Shooting", "StartingInitialConditions" -> 
                                     {x[1] == 1, y[1] == 2/100 + 91/16000}}];
Plot[{x[t], y[t]} /. sol, {t, 0, 1}, Evaluated -> True]

enter image description here

The "StartingInitialConditions" above is found by trial and error, despite some warnings generated, the solution is reliable enough:

(* Error check *)
Plot[{eqn1, eqn2} /. sol, {t, 0, 1}, Evaluated -> True]
Plot[{eqn1, eqn2} /. sol, {t, 0, 1}, PlotRange -> All]

enter image description here enter image description here


Let me add some explanations about how I found a good initial condition. As I've mentioned in the comment below, what I used is just method of exhaustion… yeah, exhaustion, sounds clumsy but it's really powerful.

First, I define the following function for the realization of the method:

SetAttributes[try, Listable]
try[i_] :=(* try[i] = *){i, Quiet@NDSolve[{eqn1 == 0, eqn2 == 0, x[0] == y[0], x[1] == 1}, 
                           {x, y}, {t, 0, 1}, Method -> {"Shooting", 
                           "StartingInitialConditions" -> {x[1] == 1, y[1] == i}}]};

If you add the try[i] = in the note into the code, Mathematica will remember all the calculated try[i] so repetitive computation can be avoided when rechecking, of course this will consume more memory and not that necessary for your case.

We don't know what value will be a proper initial condition, but it's not that hard to guess the approximate extent of it, based on some observations to the equations and boundary conditions, I guess that y[1] may be between -10 and 10. So I tried:

midpoint = 0; step = 1; 
try@Range[midpoint - 10 step, midpoint + 10 step, step] // MatrixForm

enter image description here

All of the trials fail in the middle of the domain of t, but one of them stick out to about 0.02, with y[1] == 0, which means a proper initial condition may be near 0, so I go on trying:

midpoint = 0; step = 1/10; 
try@Range[midpoint - 9 step, midpoint + 9 step, step] // MatrixForm

enter image description here

The best condition we found is still y[1] == 0 but the scope is smaller now… wait, so far I was modifying midpoint and step by hand, why not make them modified automatically?:

Clear[leftboundary, area]
(* leftboundary is used for the extraction 
   of the left boundary of the interpolating function,
   it's specifically designed for the structure of try[i]. *)
leftboundary = First@First@First@First[x /. Last@#] &;
(* area is used for the generation of possible y[1] *)
area[midpoint_, step_] := Range[midpoint - 9 step, midpoint + 9 step, step]

NestWhile[{First@SortBy[try@area[First@#, #2], leftboundary], #2/10} & @@ # &, 
          {try[0], 1/100}, leftboundary@First@# > 10^-16 &]
{513/20000, 1.24278*10^-128, 1/1000000}

OK, this time we get a even better initial condition: y[1] == 513/20000.

share|improve this answer
thanks very much – 3c. Jul 17 '13 at 13:27
by the way, do you have any idea to guess the initial condition, I tried many times, but still can not find such good initial values. thanks – 3c. Jul 17 '13 at 14:12
@user8583 What I used is just method of exhaustion… of course, this can be done automatically with the combination of Table and Nest, for this part I can add something to my answer tomorrow, but now I'd like to go to bed. :D – xzczd Jul 17 '13 at 15:15
thanks so much. I learn a lot. – 3c. Jul 22 '13 at 14:07
by the way. For your last paragraph which is talking about how to make the search automatically, I am not quite understanding what it means. such as the definition of g and the content inside nestwhile. if you have time, I appreciate your help. – 3c. Jul 22 '13 at 14:16

Your Answer


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.