# DSolve taking a very long time then goes out of memory

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]]],
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,
thetaf == 0},
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,


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.

• You should also avoid using captain letters such as R as parameter. – zhk Feb 20 '17 at 13:15
• I don't think, analytical solution is possible for your system. Do you agree? – zhk Feb 20 '17 at 13:17
• The variable in the step function should be t not x. – zhk Feb 20 '17 at 13:24
• The initial condition rd ==0 is causing 1/0. – zhk Feb 20 '17 at 13:25

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

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 == 10^(-4), thetad == Pi, thetaf == 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

• Any insight why the solution to this system of ODEs is complex? And I'm trying to find vf/vd, is there a way to do this with this approach? – OmnipotentEntity Feb 20 '17 at 15:13
• @OmnipotentEntity Apparently, because of the Sqrt's – zhk Feb 20 '17 at 15:16
• I guess that is a dumb question. 🤦 – OmnipotentEntity Feb 20 '17 at 15:17
• If you insist... – OmnipotentEntity Feb 20 '17 at 15:29
• I'm going to mark your answer as accepted and post a new question because I have fixed my code, and the question is significantly different. Thank you. – OmnipotentEntity Feb 21 '17 at 16:04