In the context of advection numerical solving, I try to implement the following recurrence formula in a time loop:
As you can see, I need the second previous time value for (j-1) and previous (j) value to compute the (j+1) time value.
I don't know how to implement this recurrence formula. Here below my attempt in Python where u
represents the array of values T
for each iteration:
l = 1
# Time loop
for i in range(1,nt+1):
# Leapfrog scheme
# Store (i-1) value for scheme formula
if (l < 2):
atemp = copy(u)
l = l+1
elif (l == 2):
btemp = copy(atemp)
l = 1
u[1:nx-1] = btemp[1:nx-1] - cfl*(u[2:nx] - u[0:nx-2])
t=t+dt
Coefficient cfl
is equal to s
.
But the results of simulation don't give fully good results. I think my way to do is not correct.
How can I implement this recurrence? i.e mostly how to store the (j-1) value in time to inject it into formula for computing (j+1) ?
In the formula:
the time index j
has to start from j=1
since we have the term T_(i,j-1)
.
So for the first iteration, we have :
T_i,2 = T_i,0 - s (T_(i+1),1 - T_(i-1),1)
Then, if In only use time loop (and not spatial loop such that way, I can't compute dudx[i]=T[i+1]-T[i-1])
, how can I compute (T_(i+1),1 - T_(i-1),1)
, I mean, without precalculating dudx[i] = T_(i+1),1 - T_(i-1),1
?
That was the trick I try to implement in my original question. The main problem is that I am imposed to use only time loop.
The code would be simpler if I could use 2D array with T[i][j]
element, i
for spatial and j
for time but I am not allowed to use 2D array in my examination.
There are few problems I see in your code. First is notation. From the numerical scheme you posted it looks like you are discretizing time with j
and space with i
using central differences in both. But in your code it looks like the time loop is written in terms of i
and this is confusing. I will use j
for space and n
for time here.
Second, this line
u[1:nx-1] = btemp[1:nx-1] - cfl*(u[2:nx] - u[0:nx-2])
is not correct since for the spatial derivatve du/dx you need to apply the central difference scheme at every spatial point of u
. Hence, u[2:nx] - u[0:nx-2]
is doing nothing like this, it is just subtracting what seems to be the solution including boundary points on the left from the solution including boundary points on the right. You need to properly calculate this spatial derivative.
Finally, the Leapfrog method which indeed takes into account the n-1
solution is usually implemented by keeping a copy of the previous time step in another variable such as u_prev
. So if you use the Leapfrog time scheme plus central difference spatial scheme, in the end you should have something like
u_prev = u_init
u = u_prev
for n in time...:
u_new = u_prev - cfl*(dudx)
u_prev = u
u = u_new
Note that u
on the LHS is to compute time n+1
, u_prev
is at time n-1
and dudx
uses u
at the current time n
. Also, you can compute dudx
with
for j in space...:
dudx[j] = u[j+1]-u[j-1]