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 have a question for solving t -> Infinity on Mathematica. First, I have a system of ODEs:

feqn1 = {x'[t] == -k1*x[t]*y[t]/h[t] + k2*z[t],
   y'[t] == 3*(-k1*x[t]*y[t]/h[t] + k2*z[t]),
   z'[t] == k1*x[t]*y[t]/h[t] - k2*z[t],
   h'[t] == -3*(-k1*x[t]*y[t]/h[t] + k2*z[t]),
   x[0] == .001, z[0] == 0, y[0] == 0.2, h[0] == .001};

Then I solved the equations by NDSolve.

NDSolve[feqn1, {x, y, z, h}, {t, 0, 1000}];

and for the second part of the problem, I have another set of ODE.

seqn1 = {x2'[t] == -k1*x2[t]*y[t]/h2[t] + k2*z[t],
         h2'[t] == -3*(-k1*x2[t]*y[t]/h2[t] + k2*z[t])}

The initial conditions are h2[0] == 1 at t = Infinity, x2 = 0.001.

Thanks.

share|improve this question
2  
You can approximate the boundary condition at infinity with a suitably large number, e.g. h2[1/$MachineEpsilon] == 1. Also, see this, among a number of other old MathGroup posts on handling such boundary conditions. –  J. M. Nov 8 '12 at 0:24
4  
An alternative is to change variables in your ODE so that infinity becomes 1 say (?) –  chris Nov 8 '12 at 8:45
    
Thanks J.M & Chris... I tried to look up the document, they used DShoot, but I couldn't find a command on Mathematica for "DShoot". And I also tried set x2[t=1000]= 0.001, but Mathematica only gave me error msg. –  DumbleKo Nov 8 '12 at 18:53
add comment

1 Answer

up vote 3 down vote accepted

The comments have made two suggestions - one that we change variables so that $[0,\infty)$ corresponds to a bounded region like $[0,1)$ and the other that we simply solve over a large interval of the form $[0,M]$. Since you've already got InterpolatingFunctions defined over $[0,1000]$ in your system, I suppose that the second suggestion probably makes the most sense. Also, as appealing as a change of variables sounds, in my experience this likely generates a singularlity at the right endpoint making anyway.

Here's your initial system; I grabbed values of $k_1$ and $k_2$ from another post.

Clear[x,y,z,h]
k1 = 0.000392;
k2 = 0.00759355499;
feqn1 = {x'[t] == -k1*x[t]*y[t]/h[t] + k2*z[t],
   y'[t] == 3*(-k1*x[t]*y[t]/h[t] + k2*z[t]),
   z'[t] == k1*x[t]*y[t]/h[t] - k2*z[t],
   h'[t] == -3*(-k1*x[t]*y[t]/h[t] + k2*z[t]),
   x[0] == .001, z[0] == 0, y[0] == 0.2, h[0] == .001};
{x,y,z,h} = {x,y,z,h} /. First[
  NDSolve[feqn1,{x,y,z,h},{t,0,1000}]
];

We'll use a shooting method to solve the boundary value problem. That is, we'll write a function shot that solves an initial value problem in terms of a parameter a, with initial conditions at the left endpoint. The idea is to choose the parameter a so that the right hand boundary condition is satisfied. Here's shot:

Clear[shot];
shot[a_?NumericQ] := Module[{},
  solution = First[NDSolve[
     {x2'[t] == -k1*x2[t]*y[t]/h2[t] + k2*z[t],
      h2'[t] == -3*(-k1*x2[t]*y[t]/h2[t] + k2*z[t]),
      h2[0] == 1, x2[0] == a},
     {h2, x2}, {t, 0, 1000}]];
  {h2Try, x2Try} = {h2, x2} /. solution;
  x2Try[1000]]
Plot[shot[a], {a, -2, 2}]

enter image description here

The plot indicates that shot is one-to-one and that there's probably a unique solution to your problem. We can find it with FindRoot:

FindRoot[shot[a] == 0.001, {a, 0}]

(* Out: {a -> -0.00479028} *)

Note that the solutions for $x_2$ and $h_2$ are already stored in the variables x2Try and h2Try. We can plot them over various intervals as below. I'm actually plotting $x_2(t)$ and $h_2(t)-1$ here since, otherwise, they look like two flat lines.

Manipulate[
 Plot[{x2Try[t], h2Try[t] - 1}, {t, 0, max}],
 {{max, 1}, 0.1, 1000}]

enter image description here

share|improve this answer
    
If need be, one could save on memory in the definition of shot[] by having the start and endpoints in the integration interval coincide: x2[1000] /. First @ NDSolve[(* equations *), {h2, x2}, {t, 1000, 1000}]. –  J. M. Nov 8 '12 at 13:31
    
@ Mark McClure: Thanks for the answer, but I'm still trying to figure out "a" because according to your answer, it seems like at t=0, x2=a, which equals to some values. But in my case, at t=0, x2=0 and t=inifinity, x2=0.001 is the boundary conditions. Maybe I misunderstood something here... –  DumbleKo Nov 8 '12 at 18:52
    
@DumbleKo At no point in your question do I see a specification for x2[0], which is a good thing! You've got a first order system in the two unknown functions h2 and x2 - correct? You've already specified that h2[0]==1. If you additionally specify that x2[0]==0, then your solution is uniquely determined so you can't expect to force the value of x2 at infinity. If you want some particular value for x2 at infinity, then you might be able to get it, if you allow flexibility at zero. That's the approach taken here. –  Mark McClure Nov 8 '12 at 20:29
    
@ Mark McClure: I see. Sorry I was confused by the my question and trying to solve this problems from different approach. But thanks for answering all my questions. They've been very helpful. –  DumbleKo Nov 8 '12 at 21:12
add comment

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.