Numerical solution of continuity equation, implicit scheme, staggered grid

[trelek2]

Hi!

I'm trying to implement an implicit scheme for the continuity equation.
The scheme is the following:
http://img28.imageshack.us/img28/3196/screenshot20111130at003.png
With \rho being the density, \alpha is a weighing constant. d is a parameter that relates the grid spacing to the velocity of flow. j is the spatial grid and n is the time grid.

The problem with this scheme is the fact that it is implicit and effectively I have no idea how to successfully implement it. I tried by doing the following:
Assume have initial time n=0 spatial (j) grid full. Also assume I know spatial boundaries (j=-1/2 and j=max) at time n+1.

Then I set j=1/2 equal to a variable x.
Next rearrange so that have j=3/2 in terms of variable x and boundary.
Do likewise for all the (half)grid points so that they can be written in terms of x. then when I reach the opposite boundary of grid (j=max) solve for x. So then I have j=1/2 and hence can substitute x into all the other equations to fill the grid points with data.

This can't work since:
When writing expression for any grid point in terms of x, i divide by 0.5d(1-\alpha). Which is roughly 0.5. So going through all the grid points I end up with 0.5^400 (400 grid points) in the denominator which is bound to kill the calculation.

Can anyone tell me how to deal with this implicit scheme?

 

[bigfooted]

Are you sure you want to work at the 1/2 nodes? Shifting everything to the left by 1/2 and you will get a scheme with unknowns at j-1,j,j+1, which will be MUCH more convenient.

First of all, assume everything at timestep n known. When n=0, you can substitute the initial conditions. Put all unknowns (the solution at n+1) to the left and all known terms to the right. This is generally a good idea for any numerical scheme.

I will now for simplicity move your scheme 1/2 node to the left. The result is:
\rho_{j} + 0.5d(1-\alpha)(\rho_{j+1}-\rho_{j-1}) = -0.5d(...)

When you have N nodes to discretize space, from j=0..N, then j=0 and j=N are on the boundaries.
Suppose j=1. We then get:
\rho_{1} + 0.5d(1-\alpha)(\rho_{2}-\rho_{0}) = -0.5d(...)

\rho_{0} is on a boundary, so if a boundary condition was imposed on x=0, then you can substitute its value (and move it to the right of the equal sign).

As long as the boundary conditions have not been substituted yet, you will have N-2 equations (j=1..N-1) for N unknowns (j=0..N).
So you have a matrix-vector system of the form [A][\rho^{n+1}]=
After substituting the boundary conditions, you should have a system of N-2 equations with N-2 unknowns.
This matrix-vector system is tri-diagonal for which very efficient algorithms exist.
If you solve this, you will get the solution at time n+1.

Hope this helps