Anatomy of a solver

Author: Bruce Wernick
Last Updated: 9 August 2011

Finding a good solution to a numerical problem is tricky. I say "a good solution", because there are always a few different approaches and in the end it's not always straight forward to decide between them.

This article illustrates the thinking behind one particular problem.

Friction is a fact of life that appears a lot in flow systems. In the early 1700, experimental measurements led to the Moody diagram that encapsulates friction in flow systems.

You can see that the friction factor depends on the dimensionless flow and the roughness. This defines a unique point on the chart. But, the friction factor relationship does not easily resolve into a single equation. In the turbulent flow region, Colebrook defined this very nice and accurate equation.

1/√f = -2 log (k /(3.7 D) + 2.51 /(Re √f)

The problem with this equation is that it can't be re-arrange to solve for f.

A good start is to try to simplify the problem by substituting for groups of variables.

Firstly, let's get rid of the 1/√f by assigning it to x. All the other constants can be grouped into C₀ and C₁.

x = -2 log₁₀ (C₀ + x C₁)

One way to proceed is to make a starting guess for f. Check the chart, f = 0.02 is a good starting point. This makes x = 7.071.

Now you can simple plug x into -2 log (C₀ + x C₁) to get the next guess. This is method is called sucessive re-substitution.

Let's try a real example. Air flowing in a duct with k/D = 0.00050 and Re=100,000.

C₀ = e/3.7 = 1.351e-4
C₁ = 2.51/Re = 2.51e-5

The first few steps look like this:

x₁ = -2*log₁₀(C₀ + x₀*C₁) = 7.0099
x₂ = -2*log₁₀(C₀ + x₁*C₁) = 7.0142
x₃ = -2*log₁₀(C₀ + x₂*C₁) = 7.0139

and that's pretty close for 3 steps.

f = 1/x² = 0.02032

Programming the solver

You can't beat Python for simplicity in solving this kind of problem.

C0 = e/3.7
C1 = 2.51/Re
x = 1/sqrt(f)
for i in range(5):
  x = -2*log10(C0 + x*C1)

That's the whole program. Notice that the code is pretty much the same as our definition of the problem. That's why I like Python so much, it's easy to experiment with alternate solutions.

With this method, we converge to within 5 decimal places in 5 steps.

But, let's look for a better solution to the problem.

We can re-arrange the starting equation by shifting all terms to the right hand side.

y = x + 2 log₁₀ (C₀ + x C₁)

So, what we want to find is the value of x when y = 0.

We refer to this as a root finding problem.

Our 1st guess x₀ = 7.071 is not the correct answer so we need to improve the guess with a bit of simple geometry.

The slope of the curve is the change in y divided by the change in x ( slope = dy / dx). Or another way round dx = dy/slope.

We use this correction to get the next guess.

x_i = x_i-1 - dx

= x_i-1 - y / slope

We need the slope of our function.

slope = 1 + 0.8686 C₁/(C₀ + x C₁)

Solving it with Python.

f = 0.02
x = 1/sqrt(f)
for i in range(5):
  y = x + 2.0*log10(C0 + x*C1)
  s = 1 + 0.8686*C1/(C0 + x*C1)
  x = x - y/s

This solution converges to 10 decimal places in 3 steps. It also converges faster with any starting guess.

The above method is called Newton-Raphson and is a well known to be a fast root finder.

Before closing the book, let's consider a potential problems.

What if the slope is zero? This means that the program will crash on the line x = x - y/s because you can't divide by zero. Even a small slope can cause problems because y/small will give a big number and therefore a big change in x. This kind of thing happens with Newton-Raphson so we need to find a way to protect against it happening.

Is there a very safe method?

Say we start with 2 guesses (x₁ and x₂), either side of the actual solution.

This means that there will be a sign change in y. The improvement of x can be found by taking the average of the last 2 guesses.

x = (x₁ + x₂) / 2

This is called the bisection method and is inherently safe. The problem is that it's much too slow.

The Python code.

f = 0.02
x = 1/sqrt(f)
x1 = x+0.1
y1 = x1 + 2.0*log10(C0 + x1*C1)
x2 = x-0.1
for i in range(50):
  x = (x1+x2)/2
  y = x + 2.0*log10(C0 + x*C1)
  if y1*y <= 0.0:
    x2 = x
  else:
    x1 = x
    y1 = y

The 1st 10 steps look like this:

7.07106781187
7.02106781187
6.99606781187
7.00856781187
7.01481781187
7.01169281187
7.01325531187
7.01403656187
7.01364593687
7.01384124937

It may be slow but does have the benefit of being safe.

By combining Newton-Raphson with bisection we can get fast convergence with safety. Essentially, we must decide if the step is bad. If so, take a bisection step.

  if badstep:
    x = (x1 + x2)/2
  else:
    x = x - y/s

The Colebrook problem has a smooth curve so we don't end up needing the bisection step.

Another problem with Newton-Raphson is that we need an equation for the slope. With a simple function like friction, we can calculate the slope but it's not always possible. Bisection doesn't need the slope but it does need 2 guesses that straddle the solution.

A much better method.

How about a method where you estimate the slope inverse once at the start. From then on, you make step-by-step corrections to the value of x and to the slope inverse.

K = 1 / slope
while 1:
  dx = -K * fo
  x = x - dx
  fx = f(x)
  if abs(fx) <= TOL: 
    return x
  dfx = fx - fo
  if abs(dfx) < TINY: 
    return x
  a = dx * K * dfx
  dK = -K * (a - dx * dx) / a
  K = K + dK
  fo = fx

This method is not much slower than Newton and only does 1 function evaluation per loop. It also starts with 1 guess and doesn't need the slope evaluation.

The 1st 4 iterations show the good progress of this method.

6.66686213905
7.0033594635
7.01395123124
7.0139611289

Broyden published this method for multi-dimensional systems. I have searched and asked a few people but not seen it simplified to 1-dimension as shown above. Maybe there is a reason that I'm not aware of but so far it has proved itself to be quite rugged.