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 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
add comment

1 Answer 1

up vote 6 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


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

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:

d = 1/10; mid = 0;
f@Range[mid - 9 d, mid + 9 d, d] // 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 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.

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
show 9 more comments

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.