Testing methods to solve Lagrangian advection¶
We'll implement a few different methods to solve the first order ODE drdt=v(r).
First of all, lets setup some useful python modules:
In [1]:
import pylab as p
Now lets define our velocity as v=(−2πy,2πx).
In [2]:
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 rn+1=rn+ΔtN∑i=0αiv(rn−i,tn−i)
where the coefficients are chosen to ensure that the method is order N, that is to say, when terms are Taylor expanded as
r(t+a)=r(t)+adrdt+a22!d2rdt2+…,
or equivalently
r(t+a)=r(t)+av(t)+a22!dvdt+…,
then the two sides match to order N.
In [3]:
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, rk=rn+ΔtN∑ianiv(rk,t+biΔt) with a final output, rn+1=rN. The explicit methods have ani=0 for i≥n, so that the levels can be found by sequential construction.
In [4]:
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 [5]:
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.
In [41]:
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:¶
Testing methods to solve Lagrangian advectionWe'll implement a few different methods to solve the first order ODE
drdt=v(r).
First of all, lets setup some useful python modules: