It was back in high school when I learned the Newton-Raphson method for computing the root of an equation using an iterative method (i.e., guessing). I belive that I learned it, applied it on whatever test came in the following weeks and set it away again – that was around 25 years ago. Last week, I got the chance to use it again, which is pretty cool, since this old method is quite clever.
My problem was that I was trying to create a little modeling widget to represent the bottom of a lake with a as a simple grid of numbers representing the “cells” elevation relative to sea level (http://sifn.bse.vt.edu/sifnwiki/index.php/Wooomm_hydroGridImpSmall). From this number I could then calculate the volume of water that would be stored in the lake at any given elevation of water surface, and conversely, if given a volume of water I could compute the surface elevation that would result. Calculating the volume at any given surface elevation is trival, you just subtract the surface value from the elevation level of the bottom of each cell and multiply the difference to the area of the cell – sum up the cells and you have your volume of water. The opposite process, back-calculating the surface elevation from a known voume, is a good bit more tricky since as the elevation changes, some cells become submereged, whereas other become ‘dried out”, and since the bottom surface is a natural land form it does not lend itself to being described simply by some sort of equation. Thus, one way to solve the surface area at a given volume is to engage in an iterative process of guessing surface elevation, and make better guesses until you hit your given volume. Or, as follows:
- make an initial guess of the surface elevation – in a realtime model you just guess whatever the value was at the last time step, but initially you need to make a guess, I just defaulted to the minimum cell bottom elevation value.
- calculate the volume of water that would be stored in the lake if the surface were at that level.
- Compare the guess to the given volume. If the guess is too high, decrease the elevation guess, if it is too low, increase the elevation guess.
This can be done fairly simply by brute force: just choosing an arbitrarily small increment for changing your elevation, and then iterating until you reach your desired volume within some error tolerance. However, this could take many steps to reach a solution, and given that thjis is to be used in simulations that will run over thousands of time steps, every small improvement in the efficiency of the solution method will amount to substantial time savings in the overall simulation. Enter the need for Newton-Raphson. Because I didn’t remember a bit about the method other than 1) it existed, and 2) it somehow enabled one to converge on a numerical solution by helping to make educated guesses, I consulted Google, who led me to Wikipedia, which had this for me: : http://en.wikipedia.org/wiki/Newton’s_method Newton’s method is a means of applying calculus for making educated guesses. The assumptions underlying this method are that the system of interest is described by a differentiable function, i.e., an equation that has a derivative. While the volume of the grid is not actually defined by a differentiable equation, the essence of a derivative is that the equations rate of change is calculable, and in this way the method to calculate volume at a given elevation is differentiable. I set up my equation as follows:
- V(h) = Σ (h – b) * a; the volume of my grid at a given surface elevation, where: h = surface elevation of lake above sea level, b = the bottom elevation of a given cell, and a = the unit area of each cell.
- f(h) = V(h) – S; this is our “differentiable” function, where S = the given storage, this function will have a slope of 0.0 (a local minima/maxima) when the elevation h yields a volume that is equal to S, the given storage.
- f'(h) = df/dh = a numerical approximation of the derivative by choosing some arbitrarily small change in h (called dh), and calculating the rate of change in volume as a function of these small differences, or numerically:
- f'(h) = df/dh = V(h-dh) – S – ( V(h+dh) – S), which simplifies to: df/dh = V(h-dh) – V(h+dh)
- f(hn+1) = hn – [ f(hn) / f'(hn) ]; the Newtons Method equation to calculate the next guess for elevation
- Applying Newtons method to obtain our next guess as follows:
- hn+1= hn – [ f(hn) / f'(h) ]
Tests saw this method yield a solution in 2-3 steps, as opposed to 50-100 using my earlier brute force attempts.