Testing methods to solve Lagrangian advection

We'll implement a few different methods to solve the first order ODE $$\tag{*} \frac{d\mathbf{r}}{dt} = \mathbf{v}\left(\mathbf{r}\right).$$ First of all, lets setup some useful python modules:

import pylab as p

Now lets define our velocity as $\mathbf{v}= (-2\pi y, 2\pi x)$.

def vel(x):
    return p.array((-2.0*p.pi*x[1], 2.0*p.pi*x[0]))

For methods we'll use the family of Adams-Bashforth explicit multistep methods and the explicit Runge-Kutta multistep methods:

Adams-Bashforth

For an Adams-Bashforth method, we have a relation $$ \mathbf{r}^{n+1} = \mathbf{r}^{n}+\Delta t \sum_{i=0}^N \alpha_i\mathbf{v}\left(\mathbf{r}^{n-i},t_{n-i}\right) $$ where the coefficients are chosen to ensure that the method is order $N$, that is to say, when terms are Taylor expanded as $$ \mathbf{r}(t+a) = \mathbf{r}(t)+a \frac{d \mathbf{r}}{dt}+\frac{a^2}{2!} \frac{d^2 \mathbf{r}}{dt^2}+\ldots, $$ or equivalently $$ \mathbf{r}(t+a) = \mathbf{r}(t)+a \mathbf{v}(t)+\frac{a^2}{2!} \frac{d \mathbf{v}}{dt}+\ldots, $$ then the two sides match to order $N$.

def ab1(x,h, fold):
    y = x + h*vel(x)
    
    return y, fold

def ab2(x,h,fold):
    f0 = vel(x)
    if len(fold)==1:
        y = x + h *(1.5*f0-0.5*fold[0])
    else:
        y = x + h*f0

    return y, [f0]

def ab3(x,h,fold):
    f0 = vel(x)
    if len(fold)<2:
        y, f0 = ab2(x,h,fold)
        if len(fold)==0:
            return y, f0
        else:
            return y, f0+fold
    else:
        f0 = vel(x)
        y = x + h *(23.0*f0-16.0*fold[0]+5.0*fold[1])/12.0
        return y, [f0]+fold[0:1]

Runge-Kutta

The multilevel Runge-Kutta methods introduce a series of sublevels, $$ \mathbf{r}_k=\mathbf{r}^{n} + \Delta t \sum_i^N a_{ni}\mathbf{v}(\mathbf{r}_k, t+b_i\Delta t) $$ with a final output, $ \mathbf{r}^{n+1}=\mathbf{r}_N $. The explicit methods have $a_{ni}=0$ for $i\geq n$, so that the levels can be found by sequential construction.

def rk2(x,h,fold):
    del fold
    k1 = vel(x)
    k2 = vel(x+h*k1)

    y = x+h*k2

    return y, []

def rk3(x,h,fold):
    del fold
    k1 = vel(x)
    k2 = vel(x+0.5*h*k1)
    k3 = vel(x+h*(2.0*k2-k1))

    y = x+h*(k1+4.0*k2+k3)/6.0

    return y, []

def rk4(x,h,fold):
    del fold
    k1 = vel(x)
    k2 = vel(x+0.5*h*k1)
    k3 = vel(x+0*5*h*k2)
    k4 = vel(x+h*k3)

    y = x+h*(k1+2*k2+2*k3+k4)/6.0

    return y, []

Plotting

We'll define a sort matplotlib routine to graph the trajectories generated by the various methods:

def plot(h, meth, data, t_max, ls=':', **kwargs):

    N = int(t_max/h)
    X = p.empty((N+1,2))

    x = p.array((1.0,0.0))
    fold = []

    for i in range(100):

        x, fold = meth(x, h, fold)

    R = p.sqrt(p.sum(x**2))
    Theta = p.arctan2(x[0],x[1])
    X[0,:] = p.array((1.0,0.0))

    for i in range(N):
        x, fold = meth(x, h, fold)
        r = p.sqrt(p.sum(x**2))/R
        theta = p.arctan2(x[0],x[1])-Theta
        X[i+1,:] = p.array((r*p.cos(theta),-r*p.sin(theta)))

    p.scatter( X[-1,0] , X[-1,1], 20, **kwargs)
    p.plot(X[:,0], X[:,1], kwargs.get('c','k')+ls)

    data[(h,meth)] = X

Finally, all that's left is to setup initial conditions and run the code.

H=[(0.01, 'o'),
   (0.02,'+'),
   (0.05,'*')]

t_max=1.0

data ={}

for h,m in H:
    fig=p.figure()
    plot(h, ab1, data, t_max=t_max, c='r',marker=m, label = 'Forward Euler')
    plot(h, ab2, data, t_max=t_max, c='b',marker=m, label = 'Adams-Bashforth 2')
    plot(h, rk2, data, t_max=t_max, c='g',marker=m, label = 'Runge-Kutta 2')
    plot(h, rk4, data, t_max=t_max, c='y',marker=m, label = 'Runge-Kutta 4')
    plot(h, ab3, data, t_max=t_max, c='m',marker=m, label = 'Adams-Bashforth 3')

    
    p.title('Timestep %4.4f'%h )
    p.axis([-1.6,1.6,-1.6,1.6])
    p.axis('equal')
    p.legend()

    t=p.linspace(0,t_max*2.0*p.pi,200)
    p.plot(p.cos(t_max*t),p.sin(t_max*t),'k-')
    p.scatter( 1.0 , 0.0, 20, c='k')  
    
    fig.savefig('test%3.2f.svg'%h)
    p.close()

Results:

Results:

Results: