VE Update: Improved Stability and Convergence

I was looking at the convergence of the first version of VE and thinking, “this can’t be right…”. When it’s coded right, Newton-Raphson typically converges extremely well, or not at all.

So looking at the code I realized I made a beginner mistake and used a one-sided numerical differentiation scheme:

sub fy
{

my $x     = shift;
my $y     = shift;
my $delta = shift;

my $f2  = f( $x, $y+$delta);
my $f1  = f( $x, $y );

return  ( $f2 - $f1 ) / $delta;
}

instead of a centred differentiation scheme:

sub fy
{
my $x     = shift;
my $y     = shift;
my $delta = shift;

my $f2  = f( $x, $y+$delta);
my $f1  = f( $x, $y-$delta );

return  ( $f2 - $f1 ) / ( 2 * $delta);
}

So I coded up the centered differencing and, presto!, the right kind of magic happened. So here, without further ado is an updated version of VE.

Version 0.2 Features

• derivatives with centered differences
• greatly improved outer loop convergence
• hooks added for Chung-like calculation when no elevation profile is given a priori
  

Convergence improvement:

The new centered differencing leads to greatly improved and much more stable convergence of the Newton outer loop.   Typically, only 3 or 4 iterations are required to converge the system to tolerances of 0.000001.

% ve.pl pt_data.csv elevation.csv 89 1.236 1.26 3.6 7.59
( crr, cda )= ( 0.005016005, 0.384299602) RESID: ( 7129890374.2886810303)
( crr, cda )= ( 0.005014958, 0.384314348) RESID: (          5.2429889725)
( crr, cda )= ( 0.005014958, 0.384314348) RESID: (          0.0000002056)

Download:

ve_0.2