DSolve taking a very long time then goes out of memory

Question

I have a rather complicated system of differential equations, and I think I'm handling it too naively. Unfortunately, I don't know enough to know how to handle this better.

Here is the code:

DSolve[{rd'[t] == 
 vd*(rd[t] - R Cos[thetad[t] - thetaf[t]])/Sqrt[
  R^2 + rd[t]^2 - 2 R rd[t] Cos[thetad[t] - thetaf[t]]], 
 thetad'[t] == 
 vd*Sqrt[(R^2 Sin[thetad[t] - thetaf[t]]^2)/(
R^2 + rd[t]^2 - 2 R rd[t] Cos[thetad[t] - thetaf[t]])]/rd[t], 
thetaf'[t] == 
vf/R*(-1)^(Re[Piecewise[{{HeavisideTheta[t], t != 0}}, 1/2]]),
rd[0] == 0,
thetad[0] == Pi,
thetaf[0] == 0},
{rd, thetad, thetaf}, t, 
Assumptions -> {t \[Element] Reals, t >= 0, R \[Element] Reals, 
 vd \[Element] Reals, vd > 0, vf \[Element] Reals, vf > 0, 
 rd[t] \[Element] Reals, rd[t] >= 0, thetaf[t] \[Element] Reals, 
 thetad[t] \[Element] Reals}]

However, I only really need the differential equation up to the first maximum of rd[t]. Because of this I have tried a version where I eliminated the *(-1)^Re[Piecewise bit, because it is irrelevant. But both versions ran for about 24 hours before crashing due to lack of memory (on a 16GB machine.)

The source of this problem is the duck and fox problem posted on 538's puzzle area a while back. Here's how I came up with this set of formulas:

I do not know the best strategy for the duck, but here's one that comes to mind. The duck will always flee the fox. And the fox will always seek to be in the closest location on the circle to the duck's position.

So essentially, we can model the location of the fox and the duck as a system of differential equations in polar coordinates with the following variables: Fox's speed v_f, duck's speed v_d, and radius of the pond R.

The fox's position can be given as:

dθ_f = v_f/R * Re{(-1)^(H(dθ_f - dθ_d))}
dr_f = 0

θ(t = 0) = 0
r(t = 0) = R

(Where H(x) is the Heaviside step function, and Re{z} denotes the real part of z.)

The duck's position is a bit more complicated:

First, we need to find the direction in which the duck is heading. For this polar coordinates are a bit of a hindrance, so we have to pop back into Cartesian for a bit.

If the duck's position is a vector given by P_d = (r_d cosθ_d * i + r_d sinθ_d * j) and the fox's is P_f = (R cosθ_f * i + R sinθ_f * j) then the angle between the two is given by:

arccos(P_d·P_f/(||P_d||||P_f||))

which simplifies down to:

arccos((r_d - R cos(θ_d - θ_f))/sqrt(R^2 + r_d^2 - 2*Rr_dcos(θ_d - θ_f)))

So

dr_d = v_d*(r_d - R cos(θ_d - θ_f))/sqrt(R^2 + r_d^2 - 2*Rr_dcos(θ_d - θ_f))
dθ_d = v_d*sqrt((R^2*sin^2(θ_d - θ_f))/(R^2+r_d^2-2*Rr_dcos(θ_d - θ_f)))/r_d
r_d(t = 0) = 0
θ_d(t = 0) = 0

Now we can solve this system of equations. And find the first maximum of rd, then set that equal to R to find the ratio vf/vd.

Answer 1

Comment

I don't think that an analytical solution is possible to your system of odes. Keeping in mind, the comments I made above, we will be able to find numerical solution.

For this I choose random values for the parameters.

R1 = 1; vd = 1; vf = 1;
Eq1 = rd'[t] ==  vd*(rd[t] - R1 Cos[thetad[t] - thetaf[t]])/
     Sqrt[R1^2 + rd[t]^2 - 2 R1 rd[t] Cos[thetad[t] - thetaf[t]]];

Eq2 = thetad'[t] == vd*Sqrt[(R1^2 Sin[thetad[t] - thetaf[t]]^2)/(R1^2 + rd[t]^2 - 
          2 R1 rd[t] Cos[thetad[t] - thetaf[t]])]/rd[t];

Eq3 = thetaf'[t] == vf/R1*(-1)^(Piecewise[{{HeavisideTheta[t], t != 0}}, 1/2]);

sys = Join[{Eq1, Eq2, Eq3, rd[0] == 10^(-4), thetad[0] == Pi, thetaf[0] == 0}];

sol = First@NDSolve[sys, {rd, thetad, thetaf}, {t, 0, 10}];

If we check the output of the solution, we will know that the solution is complex,

rd[t] /. sol /. t -> 0

0.0001 + 0. I

thetad[t] /. sol /. t -> 0

3.14159 + 0. I

thetaf[t] /. sol /. t -> 0

0. + 0. I