|
|
|
def advance_scalar(u, u_1, u_2, f, x, y, t, n, Cx2, Cy2, dt2,
V=None, step1=False):
Ix = range(0, u.shape[0]); Iy = range(0, u.shape[1])
if step1:
dt = sqrt(dt2) # save
Cx2 = 0.5*Cx2; Cy2 = 0.5*Cy2; dt2 = 0.5*dt2 # redefine
D1 = 1; D2 = 0
else:
D1 = 2; D2 = 1
for i in Ix[1:-1]:
for j in Iy[1:-1]:
u_xx = u_1[i-1,j] - 2*u_1[i,j] + u_1[i+1,j]
u_yy = u_1[i,j-1] - 2*u_1[i,j] + u_1[i,j+1]
u[i,j] = D1*u_1[i,j] - D2*u_2[i,j] + \
Cx2*u_xx + Cy2*u_yy + dt2*f(x[i], y[j], t[n])
if step1:
u[i,j] += dt*V(x[i], y[j])
# Boundary condition u=0
j = Iy[0]
for i in Ix: u[i,j] = 0
j = Iy[-1]
for i in Ix: u[i,j] = 0
i = Ix[0]
for j in Iy: u[i,j] = 0
i = Ix[-1]
for j in Iy: u[i,j] = 0
return u
D1 and D2: allow advance_scalar to be used also for \( u^1_{i,j} \):
u = advance_scalar(u, u_1, u_2, f, x, y, t,
n, 0.5*Cx2, 0.5*Cy2, 0.5*dt2, D1=1, D2=0)
|
|
|