2
$\begingroup$

I am working a problem I am working on. A boat start across a river at the point (c,0) and points its bow directly across the river to the point (0,0). The parameter a is the river current velocity, b is the boat velocity relative to the river velocity. Letting (x,y) be the position of the boat, I can show:

$$\begin{align*} \frac{dx}{dt}&=-\frac{bx}{\sqrt{x^2+y^2}}\\ \frac{dy}{dt}&=a-\frac{by}{\sqrt{x^2+y^2}} \end{align*}$$

Then I wrote:

Manipulate[
 z = NDSolveValue[{x'[t] == -b x[t]/Sqrt[x[t]^2 + y[t]^2], 
    y'[t] == a - b y[t]/Sqrt[x[t]^2 + y[t]^2], x[0] == c, 
    y[0] == 0}, {x[t], y[t]}, {t, 0, tfinal}];
 ParametricPlot[z, {t, 0, tfinal},
  PlotRange -> {{0, c}, {0, 10}},
  AspectRatio -> 1],
 {{a, 5}, 0, 10, Appearance -> "Labeled"},
 {{b, 5}, 0, 10, Appearance -> "Labeled"},
 {{c, 5}, 0, 10, Appearance -> "Labeled"},
 {{tfinal, b/c}, 0.1, 20, Appearance -> "Labeled"}]

If there was no river current, the time to cross straight over would be b/c, so I set final to b/c as a start.

However, there will be a problem as the boat approaches (0,0), as that will make the denominators above equal to zero and NDSolveValue will start to approach a problem.

Give the Manipulate a try, increasing b to see what happens. Mathematica will abort. My question: Does anyone have a suggestion how to set NDSolveValue so that the Abort problem can be avoided (or other possible suggestions so that I can explore a, b, and c and their influence on the trajectory of the boat)?

Another interesting comment. After increasing b to about 8.47, if I put the following code below my manipulate:

sol = NDSolveValue[{f'[t] == (0.4 - 0.01 s[t]) f[t], 
   s'[t] == (-0.3 + 0.005 f[t]) s[t], f[0] == 40, s[0] == 20}, {f[t], 
   s[t]}, {t, 0, 80}]

Then run this line:

Plot[sol, {t, 0, 80}]

It plots the solution, but the Manipulate program begins running again, producing this output.

NDSolveValue::mxst: Maximum number of 370917 steps reached at the point t == 0.906172238519645`. >>

Weird. Why does this happen?

$\endgroup$
5
$\begingroup$

NDSolve is smart enough to handle the denominator going to zero. However, it can't integrate past the discontinuity. Instead, let's tell it to stop when x == 0 with a WhenEvent:

DynamicModule[{z, tfinal = 0}, 
 Manipulate[
  z = NDSolveValue[{x'[t] == -b x[t]/Sqrt[x[t]^2 + y[t]^2], 
     y'[t] == a - b y[t]/Sqrt[x[t]^2 + y[t]^2], x[0] == c, y[0] == 0, 
     WhenEvent[x[t] == 0, {tfinal = t, "StopIntegration"}]}, {x[t], 
     y[t]}, {t, 0, 1000}];
  ParametricPlot[z, {t, 0, tfinal}, PlotRange -> {{0, c}, {0, 10}}, 
   AspectRatio -> 1], {{a, 5}, 0, 10, 
   Appearance -> "Labeled"}, {{b, 5}, 0, 10, 
   Appearance -> "Labeled"}, {{c, 5}, 0, 10, 
   Appearance -> "Labeled"}]]

The key bits are that we moved z and tfinal to be localized by DynamicModule (tfinal since it's no longer adjustable by Manipulate and therefore Manipulate no longer localizes it, and z just because it's good practice), and that we added a WhenEvent to set tfinal when we stop the integration:

WhenEvent[x[t] == 0, {tfinal = t, "StopIntegration"}]

This method has the benefit of not having to adjust tfinal manually anymore.

As a side note, c only changes the space scaling of the problem, and changing a and b together only changes the time scaling. You only have one (nondimensional) parameter here: the ratio b/a.

|improve this answer|||||
$\endgroup$
  • $\begingroup$ With regards to my last line, you could just fix a and c to one, and only allow adjustment of b. The minimum value for b should be 1, because otherwise, although the boat asymptotically reaches the shore, it drifts endlessly down the river! With these values, the boat always stays within the range 0 <= x <= 1 and 0 <= y <= 0.5. $\endgroup$ – 2012rcampion Mar 29 '15 at 4:28
  • $\begingroup$ Second note, if you're just interested in the path of the boat and not it's speed, you can solve for y as a function of x by reducing your set of differential equations to a single one: y'[x] == y[x]/x - Sqrt[1 + (y[x]/x)^2]/b && y[1] == 0 (which has the solution y == -x Sinh[Log[x]/b]). $\endgroup$ – 2012rcampion Mar 29 '15 at 4:35
  • 1
    $\begingroup$ The maximum distance of the boat from the straight-line trajectory (y == 0) has the cool form $\sqrt{(b-1)^{b-1}/(b+1)^{b+1}}$. Neat problem! $\endgroup$ – 2012rcampion Mar 29 '15 at 4:40
  • $\begingroup$ Great answer! I tried removing the DynamicModule part and it still worked. I am not sure why that is or is not needed? $\endgroup$ – David Mar 29 '15 at 14:32
  • $\begingroup$ @David DynamicModule just localizes the values of z and tfinal so that they don't interfere with other code that happens to use the same names. Without it, the variables are just set globally, so it still works. $\endgroup$ – 2012rcampion Mar 29 '15 at 16:00
2
$\begingroup$

Why not let Mathematica solve the system exactly?

Writing down the ODEs (scaled to a = c = 1)

eq1 = x'[t] == -  b x[t]/Sqrt[x[t]^2 + y[t]^2];
eq2 = y'[t] == 1 -  b y[t]/Sqrt[x[t]^2 + y[t]^2];

and DSolve-ing them without imposing the initial condition for the time being (to help Mathematica) gives

sol = DSolve[eq1 && eq2, {x[t], y[t]}, t]

(*
Out[55]= {{y[t] -> Sinh[(
     b C[1] - Log[
       InverseFunction[(
          b Sqrt[Cosh[C[1] - Log[#1]/b]^2 #1^2] (b + 
             Tanh[C[1] - Log[#1]/b]))/(-1 + b^2) &][-b t + C[2]]])/
     b] InverseFunction[(
       b Sqrt[Cosh[C[1] - Log[#1]/b]^2 #1^2] (b + 
          Tanh[C[1] - Log[#1]/b]))/(-1 + b^2) &][-b t + C[2]], 
  x[t] -> InverseFunction[(
      b Sqrt[Cosh[C[1] - Log[#1]/b]^2 #1^2] (b + 
         Tanh[C[1] - Log[#1]/b]))/(-1 + b^2) &][-b t + C[2]]}}
*)

This looks horrible at first sight but closer inspection shows that the expression for x[t] appears in y, so that we can write y = x Sinh[C[1] - (1/b) Log[x]]. The constant of integration is determined from y(x->1) = 0 = Sinh[C[1]] -> C[1] = 0, so that the trajectory is given by

Clear[y]

y[x_, b_] = x Sinh[-(1/b) Log[x]];

Now the time dependence is given inversely, i.e. not x[t] but t[x]. More precisely, the time as a function of x is given by inverting InverseFunction:

Clear[t]

t[x_] = -((
      Sqrt[Cosh[C[1] - Log[#1]/b]^2 #1^2] (b + Tanh[C[1] - Log[#1]/b]))/(-1 + b^2) /. C[1] -> 0) &[x]

(*
Out[70]= -((Sqrt[x^2 Cosh[Log[x]/b]^2] (b - Tanh[Log[x]/b]))/(-1 + b^2))
*)

Here we have set the arbitrary constant C[2] to 0 for simplicity.

Hence, for given x the formulas tell us y as well as t. So everything is well determined (for given boat velocity relatively to the water b).

I leave the exploration of the formulas to the reader.

|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.