using index notation for finite differences

Question

I am interested in deriving finite difference equations, which means I have to play around with indexed quantities. For example, for a 2D function $f(x,y)$, I could discretize it on a grid with grid points (i,j) and a typical finite difference model might give:

$f(i,j)-g(i,j) = a(i,j)f(i+1,j)-b(i,j)f(i,j)+c(i,j)f(i+1,j+1)$

where I know $g$,$a$,$b$ and $c$. I might want to do things like rearrange the above for $f(i+1,j+1)$, or find an expression for $f(i+n,j)$ for arbitrary n, or find coefficients of $f(i,j)$ etc. I might also want to use simultaneous equations like this.

What's the best way of representing this kind of maths so that I can perform these manipulations?

---

Edit: Here's a simple example. If I have two difference equations:

$f[i,j]=fold[i,j]+a*( g[i+1,j]-g[i-1,j] )$

and

$g[i,j]=gold[i,j]+b*( f[i,j+1]-f[i,j-1] )$

in the functions $f$ and $g$, where fold,gold,a and b are known. Then I would like to be able to solve these for $g[i,j]$ in terms of the knowns and $g$ at other indices eg. $g[i+1,j]$. To do this on paper, I use the first equation to give me expressions for $f[i,j+1]$ and $f[i,j-1]$ and then insert these into the RHS of the bottom equation. Then I would like to extract coefficients of things like $g[i+1,j+1]$ from the final expression for $g[i,j]$.

How can I do that in Mathematica?

Answer 1

These look like simple linear operations on the data. One convenient way to rpoceed is to use ListConvolve and/or ListCorrelate which convolve (or correlate) a one or two-dimensional kernal across all the data. So for instance:

ListConvolve[{{1, 0}, {-2, -1}}, {{a1, b1, c1}, {a2, b2, c2}, {a3, b3, c3}, {a4, b4, c4}}]
{{-a1 - 2 b1 + b2, -b1 - 2 c1 + c2}, {-a2 - 2 b2 + b3, -b2 - 2 c2 + c3}, 
 {-a3 - 2 b3 + b4, -b3 - 2 c3 + c4}}

applies the "kernel" {{1, 0}, {-2, -1}} to any and all elements in the data. So there is no need to deal with the individual indices of the elements. As you can see, the functions work on numbers and symbols, so can be used to setup symbolic sets of equations.

Answer 2

I think I have solved it using straight forward functions, eg. for functions f0 and f1 rather than f and g:

f0[i_, j_] = Subscript[f0, i, j]

f1[i_, j_] = Subscript[f1, i, j]

writing out "explicit" forms of f0 in terms of f1 and vice-versa:

f0ex[i_, j_] = Subscript[f0old, i, j] - a (f1[i + 1, j] - f1[i - 1, j])

f1ex[i_, j_] = Subscript[f1old, i, j] - c (Subscript[f0, i + 1, ]-Subscript[f0, i - 1, j])

use one of these explicit expressions to form the eqn

eqn1 := f1[i, j] == Subscript[f1old, i, j] - c (f0ex[i + 1, j] - f0ex[i - 1, j])

and solve for the index you want

Solve[eqn1, f1[i, j]]