Curly Trajectory

It's convenient to rewrite the system as follows:

$\dot{x}=u$

$\dot{y}=v$

$\dot{u}=-30v \sin{t}$

$\dot{v}=30u \sin{t}$

subject to the initial conditions $x(0)=y(0)=0$ , $u(0)=15$ , $v(0)=10$ .

Although we've introduced more variables, making the system larger, it is now first order, and easier to handle. Note that the sub-system in $u$ and $v$ is entirely self-contained; we can solve for $(u,v)$ first, then integrate to find $(x,y)$ .

I tried a few approaches to this.

• Numerical solution using an explicit Euler method

This felt like a well-devised trap that I fell into. The method fails, terribly, even for very small step sizes.

[Edit: turns out, this wasn't a trap at all and should have worked - see the comments below]

• Numerical solution using an implicit Euler method

This worked better, though still needed small step sizes to work. The scheme is as follows:

$\mathbf{r_{k+1}}=(\mathbf{I}-h\mathbf{A_{k+1}})^{-1} \mathbf{r_k}$

where $\mathbf{r_k}=\begin{pmatrix} u_k\\ v_k \end{pmatrix}$ , $\mathbf{A_k}=\begin{pmatrix} 0 & -30\sin{t_k}\\ 30\sin{t_k} & 0 \end{pmatrix}$ , $t_k = kh$ , and $h$ is the step-size. This is not as bad to implement as I made it look!

Once the $u$ and $v$ values were calculated up to $t=12$ , I just used numerical integration to get $x$ and $y$ (the position of the particle after this time).

The results for a few $h$ -values are below:

These values seem to be converging, but slowly. Quite a handy trick here is to extrapolate from these data points using the ratio of the differences in the successive answers. Using this gave an estimate of $16.97875\ldots$ for the answer.

• Partial analytic solution

And this is what I should have done first. Going back to the system of equations, we have

$\dot{u}=-30v \sin{t}$

$\dot{v}=30u \sin{t}$

Dividing through, $\frac{\dot{u}}{\dot{v}}=-\frac{v}{u}$

Simplifying, $\dot{u}u+\dot{v}v=0$

so that $u^2+v^2=\text{const.}$

From the initial values, the constant is $325$ .

We can now substitute $v=\sqrt{325-u^2}$ to give

$\dot{u}=-30\sqrt{325-u^2} \sin{t}$ , which is a separable equation. With the initial condition $u(0)=15$ , the solution is

$u(t)=-5\sqrt{13} \sin {\left( -30 \cos(t)+30-\sin^{-1} {\left(\frac{3}{\sqrt{13}}\right)} \right)}$

Likewise, $v(t)=5\sqrt{13} \sin {\left( -30 \cos(t)+30+\sin^{-1} {\left(\frac{2}{\sqrt{13}}\right)} \right)}$

Integrating these numerically, we find that $x(12)=-13.0505\ldots$ and $y(12)=10.861\ldots$ , giving the particle's distance from the origin as $\boxed{16.9787\ldots}$ .

Simulation code (Python). This uses explicit Euler integration with a very small time step.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

import math

dt = 10.0**(-7.0)
t = 0.0

m = 1.0

x = 0.0
y = 0.0

xd = 15.0
yd = 10.0

xdd = 0.0
ydd = 0.0

while t <= 12.0:

    x = x + xd * dt
    y = y + yd * dt

    xd = xd + xdd * dt
    yd = yd + ydd * dt

    Fx = -30.0 * yd * math.sin(t)
    Fy = 30.0 * xd * math.sin(t)

    xdd = Fx/m
    ydd = Fy/m

    t = t + dt

print math.hypot(x,y)  # Returns 16.98

I do not have any solution to add to this problem. Even I solved it numerically. I was just playing around with the initial conditions and I was surprised to see the genesis of some very beautiful looking solutions. I have attached a screenshot.

What I find really remarkable is that solving a simple set of equations can yield such patterns. There is indeed beauty hidden in mathematics. Maybe this tells us more than what we can appreciate with our current state of knowledge. Good examples of such patterns are the Lorenz strange attractor and the Mandelbrot set.

As always, a very nice problem!

Here is my code, which uses Explicit Euler Integration. Given the nature of the differential equations that describe the system, I don't think Explicit Euler is the best method. The code is actually pretty inaccurate for even slightly greater timesteps; the timestep which will yield a satisfactory result is $10^{-6}$ . It needs hardware acceleration; CUDA python is the best way to go. At the end of the code, the results are shown at various values of deltaT, as are the times taken to compile the code at the values of deltaT.

The code takes a long time to compile as the step size increases. The error in the result is dependent on step size. The Euler method is inaccurate and the convergence of the result is very slow. It is also very time consuming to compute such small values of deltaT on the CPU, and get an inaccurate result nevertheless. A better method which has $O(n^4)$ error would be Runge-Kutta.

Edit: I have tried it with a modified Euler technique below as well:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

import math

x = 0
y = 0


time = 0
deltaT = 10**-6


xDot = 15
yDot = 10




while time <= 12:
    xDotDot = -30*yDot*math.sin(time)
    yDotDot = 30*xDot*math.sin(time)

    xDot = xDot + xDotDot*deltaT
    yDot = yDot + yDotDot*deltaT

    x = x + xDot*deltaT 
    y = y + yDot*deltaT 


    time = time + deltaT


print("x: " +str(x))
print("y: " +str(y))
print("distance: " +str(math.hypot(x,y)))

#Results: 

#deltaT = 10**-1: distance = 6.100549268950116e+39 (error is dependent on size of deltaT; 0.1 too large a step size for Euler)
#deltaT = 10**-2: distance = 672596301067.0984 (still very large)
#deltaT = 10**-3: distance = 71.69830313999711 (getting into 2 digits)
#deltaT = 10**-4: distance = 19.218083697766765 (slightly more accurate)
#deltaT = 10**-5: distance = 17.186272731008557 (reasonably accurate)
#deltaT = 10**-6: distance = 16.999350040357143 (a good result)
#deltaT = 10**-7: distance = 16.98(approx) (results converge towards 16.98)


#specs: Intel i7-9750H laptop processor; 16GB 3200 Mhz RAM
#time taken to compile code for values of deltaT

# <1 sec - 10**-1
# <1 sec - 10**-2
# <1 sec - 10**-3
# <1 sec - 10**-4
#  2 sec - 10**-5
#  14 sec - 10**-6
#  2 minutes - 10**-7