# 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[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?

-

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]


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[f]
SetAttributes[f, Listable]
f[i_] :=(* f[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 f[i] = in the note into the code, Mathematica will remember all the calculated f[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:

d = 1; mid = 0;
f@Range[mid - 10 d, mid + 10 d, d] // 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[1] == 0, which means a proper initial condition may be near 0, so I go on trying:

d = 1/10; mid = 0;
f@Range[mid - 9 d, mid + 9 d, d] // MatrixForm


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

Clear[g, h]
(* g is used for the extraction of the left boundary of the interpolating function,
it's designed for the output of f[i]. *)
g = First@First@First@First[x /. Last@#] &;
(* h is used for the generation of possible y[1] *)
h[midpoint_, interval_] := Range[midpoint - 9 interval, midpoint + 9 interval, interval]

NestWhile[{First@SortBy[f@h[First@#[[1]], #[[2]]], g], #[[2]]/10} &,
{f[0], 1/100}, g@#[[1]] > 10^-16 &]


{513/20000, 1.24278*10^-128, 1/1000000}

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

-
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