# 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:

In :
import pylab as p


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

In :
def vel(x):
return p.array((-2.0*p.pi*x, 2.0*p.pi*x))


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

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$.

In :
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)
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+5.0*fold)/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.

In :
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:

In :
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,x)
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,x)-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.

In :
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:¶ 