6
$\begingroup$

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?

Improve this question
$\endgroup$
0

1 Answer 1

Reset to default
13
$\begingroup$

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[try]
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.

Improve this answer
$\endgroup$
11
  • 1
    $\begingroup$ 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 $\endgroup$
    – 3c.
    Jul 17, 2013 at 14:12
  • 1
    $\begingroup$ @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 $\endgroup$
    – xzczd
    Jul 17, 2013 at 15:15
  • $\begingroup$ 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. $\endgroup$
    – 3c.
    Jul 22, 2013 at 14:16
  • $\begingroup$ @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… $\endgroup$
    – xzczd
    Jul 22, 2013 at 15:14
  • $\begingroup$ 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 $\endgroup$
    – 3c.
    Oct 18, 2013 at 18:30

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.