# using index notation for finite differences

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?

• Is this question about the software Mathematica or about math in general? Commented Oct 22, 2013 at 17:54
• Nasser has quite a few demos and functions that do this kind of stuff... here's an example from chat.
– rm -rf
Commented Oct 22, 2013 at 17:55
• It's about Mathematica. I can do it on paper but it requires pages and pages. Commented Oct 22, 2013 at 17:57
• Have you read the extensive NDSolve documentation? In particular tutorial/NDSolvePDE and stuff about method of lines. Commented Oct 22, 2013 at 20:50
• I haven't read the entire NDSolve doc, but from my experience that is all about solving PDE's numerically. That's not what I want to do, I want to do numerical analysis, which involves manipulating the difference equations themselves (rather than solving difference equations). Commented Oct 23, 2013 at 13:34

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.

• Thanks for your answer but I can't see how to apply these to the problem. I've edited my question with a simple example, could you show how that is possible with ListConvolve? It may be that straight forward functions will work but I haven't managed to do this yet. Commented Oct 23, 2013 at 13:32

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]]