Numerically solving Differential Equation using Störmer-Verlet

40 Views Asked by At

For uni I've got this assignment where I need to numerically solve the differential equation: $$\begin{equation} \begin{cases} x' = y \\ y' = -x \end{cases} \end{equation}$$ with $x_0 = 1, y_0 = 0$ numerically. Specifically I need to use the Störmer-Verlet method of solving this. To do this I first needed to rewrite this DE into a second degree DE in the format $\ddot{x} = f(x)$. I did this using a method described in my book, and got $\ddot{x} = -x$. Now, using Störmer-Verlet i get the following description of $x_n$: $$x_n = (2-(dt)^2)x_{n-1} - x_{n-2}$$

This is all fine and dandy, but for my plot I still need a discription of $y_{n}$. My first thought was that I could use $$y_n = \frac{x_n - x_{n-1}}{dt}$$ as according to the original DE, $y = x'$ was still true. However, when implementing the following code I get the plot shown below (it's the purple plot labeled SV).

def störmer_verlet(x_0, x_1, y_0, t_delta, n_iterations):
    y_1 = (x_1-x_0)/t_delta
    result_x = [x_0, x_1]
    result_y = [y_0, y_1]
    x_i = x_1
    y_i = y_1
    for i in range(2, n_iterations+1):
        x_imin1 = result_x[-1]
        x_imin2 = result_x[-2]
        x_i = (2-t_delta) * x_imin1 - x_imin2
        y_i = (x_i - x_imin1)/t_delta
        result_x.append(x_i)
        result_y.append(y_i)
    return result_x, result_y

t_delta = 0.1
x_0 = 1
y_0 = 0
x_1 = x_0 + y_0*t_delta + (t_delta**2)/2 * (-x_0)
SV_x, SV_y = störmer_verlet(x_0, x_1, y_0, t_delta, n)

fig, ax = plt.subplots()
fig.set_size_inches(7, 7)
ax.set_aspect('equal')
f = lambda t: (cos(t), sin(t))
t_values = np.linspace(0, 6.3, 64)
x_real = list(map(lambda t: cos(t), t_values))
y_real = list(map(lambda t: sin(t), t_values))
ax.plot(x_real, y_real, label='actual')
ax.plot(SV_x, SV_y, label='SV', linestyle='dashdot')
" OTHER METHODS "
ax.legend()

This is obviously not as accurate as this Störmer-Verlet method is supposed to be. Its obvious enough that the error is present mostly in the y direction. Could anyone help me towards another way of generating $y_n$ which is more accurate? Thanks in advance!!

enter image description here