# How can I deal with a non-numerical value for a derivative at $t = 0$ when using NDSolve?

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 == y, x == 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?

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 == y, x == 1}, {x, y}, {t, 0, 1},
Method -> {"Shooting", "StartingInitialConditions" ->
{x == 1, y == 2/100 + 91/16000}}];
Plot[{x[t], y[t]} /. sol, {t, 0, 1}, Evaluated -> True] 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]  Edit

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:

ClearAll[try]
SetAttributes[try, Listable]
try[i_] :=(* try[i] = *){i, Quiet@NDSolve[{eqn1 == 0, eqn2 == 0, x == y, x == 1},
{x, y}, {t, 0, 1}, Method -> {"Shooting",
"StartingInitialConditions" -> {x == 1, y == 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 may be between -10 and 10. So I tried:

midpoint = 0; step = 1;
try@Range[midpoint - 10 step, midpoint + 10 step, step] // MatrixForm 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 == 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 The best condition we found is still y == 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 *)
area[midpoint_, step_] := Range[midpoint - 9 step, midpoint + 9 step, step]

NestWhile[{First@SortBy[try@area[First@#, #2], leftboundary], #2/10} & @@ # &,
{try, 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 == 513/20000.

• 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, 2013 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 Jul 17, 2013 at 15:15
• 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, 2013 at 14:16
• @user8583 Oh, that's a combination of "pure function" and "list manipulation", which are two of the most important parts of the core language of Mathematica, you can have a look at the documents of Function and those about list manipulation. Also, you can read through this post. If you still feel confused after check all these pages, I don't mind to add more explanations to my answer… Jul 22, 2013 at 15:14
• I am trying to understand your last paragraph about how to search the initial condition automatically, would you mind dividing it into a simpler way? sorry for that, thanks
– 3c.
Oct 18, 2013 at 18:30