# Fighting the Current

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?

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.

• 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. – 2012rcampion Mar 29 '15 at 4:28
• 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]). – 2012rcampion Mar 29 '15 at 4:35
• 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! – 2012rcampion Mar 29 '15 at 4:40
• Great answer! I tried removing the DynamicModule part and it still worked. I am not sure why that is or is not needed? – David Mar 29 '15 at 14:32
• @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. – 2012rcampion Mar 29 '15 at 16:00

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.