The document is an introduction to the newton raphson solver that currently resides in TIME.

Why was it built?

The Wetlands module rendered a cluster of wetland nodes and conveyance links into a series of mass balance equations (the exact equations can be found in the wetlands SRG entry). It then tries to find the zero (or at least minimum) root of those equations - the point at which everything works.

The first draft of the Wetlands function used a very simple root finder. At the time there was not an implementation for the Multivariate Newton Raphson solver that had an appropriate licence.

In 2015 the Wetlands module was revisited, and it was found that Math.Net.Numerics had an appropriate licence and a Univariate Newton Raphson Solver. Math.Net.Numerics also had a lot of the underlying Matrix functions needed to make the Multivariate solver work. It was decided to build one.

How does it work?

How a Univariate Newton Raphson Solver Works

A newton Raphson solver works (conceptually) by looking at the local data.

Say our formula was F=x-4
If we set x=3, we can see that the value of F2= -1. We can also see that Fx=1, ie the gradient of the line is 1. For every 1 we add to x, we add one to F.

we can see how far we need to go to get to the zero by taking our current value and dividing it by the gradient. FxFx =-11 =-1

We then subtract this from our starting point of x=3, and get x=3--1=4

Unknown Functions

Say we didn’t know the function F(x), but can still evaluate it.

We guess a value - let’s say x=3. After evaluation we discover that F(3)=31. While we cannot differentiate an unknown function, we can find a nearby point and take the gradient between them. F(3.1)=33.791. That means there is a gradient of 2.791/0.1 = 27.91.

we divide the value of F(3)=31by the gradient and get 31/27.91=1.11.

We subtract this from our original guess and get x=1.89

That’s not zero, so we keep going. We use this as our new guess and repeat

at x - 1.89, F=010.74, F’=11.87, m=11.29, and new x =0.94. Again, this is not zero, so repeat. This is shown in the table below.

x	F	F'	m	subtract	newx
3.00	31.00	33.79	27.91	1.11	1.89
1.89	10.74	11.87	11.29	0.95	0.94
0.94	4.82	5.12	2.93	1.65	-0.71
-0.71	3.64	3.77	1.31	2.77	-3.48
-3.48	-38.30	-34.76	35.38	-1.08	-2.40
-2.40	-9.86	-8.20	16.60	-0.59	-1.81
-1.81	-1.91	-0.98	9.27	-0.21	-1.60
-1.60	-0.11	0.61	7.23	-0.02	-1.59
-1.59	0.01	0.71	7.09	0.00	-1.59
-1.59	0.00	0.71	7.09	0.00	-1.59

The function was actually F=x3+4, and the true root is at 3-4 =-1.5874

If we knew the function we could simply have said Fx=3x2. This would have given us the following.

x	F	m	subtract	newx
3.00	31.00	27.00	1.15	1.85
1.85	10.35	10.29	1.01	0.85
0.85	4.60	2.15	2.15	-1.30
-1.30	1.80	5.07	0.36	-1.66
-1.66	-0.54	8.22	-0.07	-1.59
-1.59	-0.02	7.59	0.00	-1.59
-1.59	0.00	7.56	0.00	-1.59

It can be seen that this follows a similar path, but much more quickly.

An interesting case happens if our estimate of m is too coarse. If we used F’=F(x+1) - ie incremented by 1, we would get the following

x	F	F'	m	subtract	newx
3.00	31.00	68.00	37.00	0.84	2.16
2.16	14.11	35.62	21.51	0.66	1.51
1.51	7.42	19.74	12.33	0.60	0.90
0.90	4.74	10.91	6.17	0.77	0.14
0.14	4.00	5.47	1.46	2.73	-2.60
-2.60	-13.55	-0.09	13.46	-1.01	-1.59
-1.59	-0.04	3.79	3.83	-0.01	-1.58
-1.58	0.04	3.80	3.77	0.01	-1.59
-1.59	-0.04	3.79	3.83	-0.01	-1.58
-1.58	0.04	3.80	3.77	0.01	-1.59

In this case the solution cannot converge. The estimate is simply too coarse to find the exact number

No Zero Root

So, what about the function F=x2+4?

This function has a minimum point of 4.

x	F	m	subtract	newx
3.00	13.00	6.00	2.17	0.83
0.83	4.69	1.67	2.82	-1.98
-1.98	7.93	-3.97	-2.00	0.02
0.02	4.00	0.03	119.51	-119.49
-119.49	14,281.75	-238.98	-59.76	-59.73
-59.73	3,571.44	-119.46	-29.90	-29.83
-29.83	893.86	-59.66	-14.98	-14.85
-14.85	224.47	-29.70	-7.56	-7.29
-7.29	57.14	-14.58	-3.92	-3.37
-3.37	15.36	-6.74	-2.28	-1.09
-1.09	5.19	-2.18	-2.38	1.29
1.29	5.65	2.57	2.20	-0.91
-0.91	4.83	-1.82	-2.65	1.74
1.74	7.02	3.47	2.02	-0.28
-0.28	4.08	-0.57	-7.20	6.91
6.91	51.81	13.83	3.75	3.17

The algorithm closes in on zero. As it reaches the zero point, the gradient moves towards 0 as well (a constant). This causes the next guess to be wildly off. The algorithm is very unlikely to converge.

We still get some accurate guesses - our best is x=0.02at which point F=4.

For the purposes of the algorithm, functions with constant areas (m approaches 0) are bad. The algorithm has little understanding of how to improve the answer.

Multivariate Newton Raphson Solver

Solving for multiple variables is just a variation on the Univariate Newton Raphson Solver.

New example:

There are three functions

F0=3x0+2x1

F1=2x0-x22

F2=x2-4

Note there is no relationship between the subscripts for F and x ( F0need not be related to x0, for instance).

The three functions then become the Matrix F.

The three variables become x.

If we have a set of values for x, V, then we can calculate FV.

so if x0=1 and x1=2 and x2=3 we can calculate the value of F at these points

F0=3(1)+2(2)=7

F1=2(1)-(3)2=-7

F2=(3)-4=-1

These are nowhere near zero.

We can use a Multivariate Newton Raphson solver to find what we should do to get closer to zero.

The steps used in the Univariate Newton Raphson Solver were:

Have a guess V
Find the value of F at V
Find the gradient of F at V
New V = V-(value of F at V/gradient of F at V)
Repeat until converged

We use the same steps here. The differences are that:

F is a matrix.
The value of F(V) is a vector
The gradient of F(V) is done by finding the Jacobian
The new value of V is a vector.

Though conceptually the same, this does have a few extra constraints.

Underspecified and Overspecified

In order to solve a system of simultaneous equations using a matrix, we must have as many variables as we do functions. Not more, not less. If we have more functions it is overspecified. If it has less functions it is underspecified.

If any function is a constant at V, then at V this will cause the problem to fail. This makes sense even in a univariate case: if F(x) returns the same value when we modify x, then it does not help us to solve a problem.

Imagine the following simple example

F0^=2x0+3x1

F1=x1+1

We can then say

x1=F1- 1

We’re looking for the root, which means F1=0

x1=-1

Which means

F0^=2x0+3(F1-1)

0=2x0+3(-1)

0=2x0-3

3=2x0

x0 =1.5

However, if we have

F1=1

F1will always be equal to 1. It gives us no clues as to how to solve Fo

However, if we have a constant function we can remove AND a variable we can remove, we can still end up with a solvable problem

F0^=2x0+3x1

F1=x1+1

F2=0x2+1

Here, F2is a constant. It does nothing to help us solve x2but on the other hand, nothing uses x2

We can simply remove it.

F0^=2x0+3x1

F1=x1+1

Optimisations

There are some things we can do to optimise the solver. These are a little dirty, but help.

Check if we are already have the zero root at the start. If we do, stop.
If we are oscillating around the zero root, then the method is too coarse. Dampen the oscillations by just finding the midpoint between the best results we have.

When Wetlands uses the Newton Raphson Solver

When wetlands builds a cluster of cells (see the wetland SRG) it creates a series of functions representing the mass balance of each cell.

Imagine we have cells A, B and C.

there is a conveyance link between A and B, and between B and C.

Each cell has a height, a, b and c.

There are 3 single value functions for mass balance - one for each cell

F0(a), F1(b), and F2(c)

There are also functions for the conveyance between each cell

F3(a,a), F4(a,b), F5(a,c), F6(b,a), F7(b,b), F8(b,c), F9(c,a), F10(c,b), F11(c,c)

These give the three full cell functions

FA=F0(a)+F3(a,a)+F4(a,b)+F5(a,c)

FB=F1(b)+F6(b,a)+F7(b,b)+F8(b,c)

FC=F2(c)+F9(c,a)+F10(c,b)+F11(c,c)

Obviously there is no flow between a node and itself, so we can remove all those function components

FA=F0(a)+F4(a,b)+F5(a,c)

FB=F1(b)+F6(b,a)+F8(b,c)

FC=F2(c)+F9(c,a)+F10(c,b)

In this case there is no conveyance link between A and C, meaning F5(a,c) and F9(c,a) will always be zero and can be removed

FA=F0(a)+F4(a,b)

FB=F1(b)+F6(b,a)+F8(b,c)

FC=F2(c)+F10(c,b)

Flows are symmetrical - F(n,m)=-F(m,n)

FA=F0(a)+F4(a,b)

FB=F1(b)-F4(a,b)+F8(b,c)

FC=F2(c)+F8(b,c)

Finally, in this case cell A is a Wetlands Hydraulic Connector. One of their properties is that they always mass balance. This means F0(a)=0

FA(a,b)=F4(a,b)

FB(a,b,c)=F1(b)-F4(a,b)+F8(b,c)

FC(c)=F2(c)+F8(b,c)

We now have the functions to solve using the Newton Raphson Multivariate Solver

What about if the Hydraulic Connector’s Conveyance Link doesn’t pass water.

Browser not supported