Solve system of ordinary differential equations that doesn't have an initial condition (t=0), but has an inifinity condition (t=infinity)?

Question

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.

Answer 1

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}]

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}]