Let's do a little (just a little) algebra.
Let x be the value in the central square; let x_i, i = 1, .., 8 index the values in the neighboring squares; and let r be the topographic ruggedness index. This recipe says r^2 equals the sum of (x_i - x)^2. Two things we can compute easily are (i) the sum of the values in the neighborhood, equal to s = Sum{ x_i } + x; and (ii) the sum of squares of the values, equal to t = Sum{ x_i^2 } + x^2. (These are focal statistics for the original grid and for its square.)
Expanding the squares gives
r^2 = Sum{ (x_i - x)^2 }
= Sum{ x_i^2 + x^2 - 2*x*x_i }
= Sum{ x_i^2 } + 8*x^2 - 2*x*Sum{x_i}
= [Sum{ x_i^2 } + x^2] + 7*x^2 - 2*x*[Sum{ x_i } + x - x]
= t + 7*x^2 - 2*x*[Sum{ x_i } + x] + 2*x^2
= t + 9*x^2 - 2*x*s.
For example, consider a neighborhood
1 2 3
4 5 6
7 8 9
Here, x = 5, s = 1+2+...+9 = 45, and t = 1+4+9+...+81 = 285. Then
(1-5)^2 + (2-5)^2 + ... + (9-5)^2 = 16 + 9 + 4 + 1 + 1 + 4 + 9 + 16 = 60 = r^2
and the algebraic equivalence says
60 = r^2 = 285 + 9*5^2 -2*5*45 = 285 + 225 - 450 = 60, which checks.
The workflow therefore is:
Given a DEM.
Compute s = Focal sum (over 3 x 3 square neighborhoods) of [DEM].
Compute DEM2 = [DEM]*[DEM].
Compute t = Focal sum (over 3 x 3 square neighborhoods) of [DEM2].
Compute r2 = [t] + 9*[DEM2] - 2*[DEM]*[s].
Return r = Sqrt([r2]).
This consists of 9 grid operations in toto, all of which are fast. They are readily carried out in the raster calculator (ArcGIS 9.3 and earlier), the command line (all versions), and Model Builder (all versions).
BTW, this is not an "average elevation change" (because elevation changes can be positive and negative): it is a root mean square elevation change. It is not equal to the "topographic position index" described at http://arcscripts.esri.com/details.asp?dbid=14156 , which (according to the documentation) equals x - (s - x)/8. In the example above, the TPI equals 5 - (45-5)/8 = 0 whereas the TRI, as we saw, is Sqrt(60).
It seems that all the values in slopes are rounded to the nearest 0,5. So 0,401 m / 2.8284 m = 0,14177, would instead be 0,5/3 = 16,6667. Or if distance is 2 metres, then 0,5/2 = 25 %. In the manual it says the values are rounded in FLAT areas, which seems to be a little off considering it does that where slope is 25 %. Another remark is that in flat areas (lakes) the drop percentages are indeed continuous!
Best Answer
I know this is not addressing your query directly, i.e., I can't give you info on calculating slope yourself. But you can calculate slope using a built in tool, well, a plugin.
This link: http://www.gistutor.com/quantum-gis/20-intermediate-quantum-gis-tutorials/48-quantum-gis-qgis-raster-based-terrain-analysis-techniques.html shows you how.
If you want to calculate it yourself then let me know and we can find out how to do raster calculations manually. Here's how to calculate it: http://en.wikipedia.org/wiki/Geographic_information_system#Slope_and_aspect