import numpy as np
import scipy as sc
import scipy.sparse.linalg
from numpy.linalg import solve
import matplotlib.pyplot as plt

central_difference=False
# set simulation parameters
h=0.25
L=1.0
n = int(round(L/h))
Tb=25 #rhs
sigma=100
k=1.65 
beta = sigma*L**2/k

y = np.arange(n+1)*h

def analytical(x):
    return beta*(1-x*x)/2+Tb
def tri_diag(a, b, c, k1=-1, k2=0, k3=1):
    """ a,b,c diagonal terms
        default k-values for 4x4 matrix:
        | b0 c0 0  0 |
        | a0 b1 c1 0 |
        | 0  a1 b2 c2|
        | 0  0  a2 b3|
    """
    return np.diag(a, k1) + np.diag(b, k2) + np.diag(c, k3)
# defina a, b and c vector
a=np.ones(n-1)
b=..
c=..

if central_difference:
    c[0]= ...
else:
    b[0]=...

A=tri_diag(a,b,c)
print(A) # view matrix - compare with N=4 to make sure no bugs
# define rhs vector
d=...
#rhs boundary condition
d[-1]=...

Tn=np.linalg.solve(A,d)
print(Tn)

The correct solution for \( L=1 \) m, and \( h=1/4 \), is: $[T_0,T_1.T_2,T_3]$=[55.3030303 , 53.40909091, 47.72727273, 38.25757576] (central difference) and $[T_0,T_1.T_2,T_3]$=[62.87878788, 59.09090909, 51.51515152, 40.15151515] (forward difference)

Exercise 4.3: Solve the full heat equation

Part 1. Replace the time derivative in equation \eqref{eq:nlin:heateq} with

and show that by using an implicit formulation (i.e. that the second derivative with respect to \( x \) is to be evaluated at \( T(t+\Delta t)\equiv T^{n+1} \)) that equation \eqref{eq:nlin:heateq} can be written

where \( \alpha\equiv\rho c_p/k \).

Part 2.

Use the central difference formulation for the boundary condition and show that for four nodes we can formulate equation \eqref{eq:lin:imp} as the following matrix equation

Part 3. Assume that the initial temperature in the concrete is $25^\circ$C, $\rho$=2400 kg/m$^3$, a specific heat capacity \( c_p= \) 1000 W/kg K, and a time step of \( \Delta t=86400 \) s (1 day). Solve equation \eqref{eq:lin:heatfull}, plot the result each day and compare the result after 50 days with the steady state solution in equation \eqref{eq:nlin:heatsol}.

Exercise 4.4: Using sparse matrices in python

In this part we are going to create a sparse matrix in python and use scipy.sparse.linalg.spsolve to solve it. The matrix is created using scipy.sparse.spdiags.

Part 1. Extend the code you developed in the last exercises to also be able to use sparse matrices, by e.g. a logical switch. Sparse matrices may be defined as follows

import scipy.sparse.linalg

#right hand side
# rhs vector
d=np.repeat(-h*h*beta,n)
#rhs - constant temperature
Tb=25
d[-1]=d[-1]-Tb
#Set up sparse matrix
diagonals=np.zeros((3,n))
diagonals[0,:]= 1
diagonals[1,:]= -2  
diagonals[2,:]= 1
#No flux boundary condition
diagonals[2,1]= 2
A_sparse = sc.sparse.spdiags(diagonals, [-1,0,1], n, n,format='csc')
# to view matrix - do this and check that it is correct!
print(A_sparse.todense())
# solve matrix
Tb = sc.sparse.linalg.spsolve(A_sparse,d)

# if you like you can use timeit to check the efficiency
# %timeit sc.sparse.linalg.spsolve( ... )

• Compare the sparse solver with the standard Numpy solver using %timeit, how large must the linear system be before an improvement in speed is seen?

CO$_2$ diffusion into aquifers

The transport of CO$_2$ into aquifers can be described according to the diffusion equation

where \( C(z,t) \) is the concentration of \co\ as a function of depth (\( z \)) and time \( t \), and \( K(z) \) is the diffusion constant of \co\ as a function of depth. This equation can be discretized using standard techniques, to help in that respect consider figure 21.

In the following we will assume that there are only four nodes (\( i=0\ldots 3 \)) in the physical domain, and two ghost nodes \( i=-1 \), and \( i=4 \). There are many ways to attack this problem, but in the following we will borrow ideas from Finite Volume. Finite volume methods is a way of discretizing equations such that we conserve mass. The diffusion equation as it is derived in figure 19, express that the flux of something (heat, particles, etc) leaving the box surface minus the flux entering the surface of the box is equal to the rate of change of something inside the box. We can formulate this mathematically as:

The notation \( i\pm1/2 \), means that the flux is to be evaluated at the surface of the box (i.e. halfway between the red dots in figure 21). \( K(z) \) is the diffusion constant, and it is known everywhere, so this is simple to evaluate at the surface. The concentrations are only known at the center of each box, the red dots in figure 21. The derivative of the concentration can be evaluated using the central difference formula (remember that the distance between the red dot and edge of the box is \( h/2 \)), hence

notice that we have discretized the time derivative, and that we have introduced \( n \) to indicate the time step. On the right hand side there are is no time indicated, it turns out that we have a choice to put time step \( n \) or \( n+1 \) on the concentrations on the right hand side. If we put \( n \) the scheme is said to be explicit, if we put \( n+1 \), the scheme is implicit. Implicit schemes are stable compared to explicit schemes, whereas explicit schemes has slightly higher numerical accuracy (TO DO 1: show this!) . In general we can write

hence if \( \theta=1 \) the scheme is explicit, if \( \theta=0 \) the scheme is implicit, and if \( \theta=1/2 \), the scheme is called the Crank-Nicolson method. The first and last boundary are special, let us first consider the \( i=0 \), this is where the sea is in contact with the \co\ in the atmosphere, and the flux is \( k_w(C_0-C_{eq}) \), hence

For the last block the flux is zero towards the seafloor, and equation \eqref{eq:lin:co2fv23} can be written

For the blocks \( i=1\ldots2 \), we can collect all terms with \( n+1 \) on one side and terms with \( n \) on the other side and rewrite equation \eqref{eq:lin:co2fv23}

where \( \alpha\equiv\Delta t/h^2 \). Next, we want to write down the corresponding matrix equations for four grid nodes as indicated in figure 21. Notice that we need to use the equations in figure 21, for \( C_{-1} \), and \( C_4 \). The left and right hand coefficient matrix \( \mathbf{L} \), and \( \mathbf{R} \) are given as

respectively. Introducing \( \mathbf{S}=\left[k_wC_{eq}\Delta t/h,0,0,0\right]^T \), we can finally write the diffusion equation \eqref{eq:lin:co2_diff} as

More stuff to do:

1. Assume zero flux over the air water interface ($k_w$=0), show from the equations above that if we start with a uniform concentration in the sea ($\mathbf{C}^n$=constant) that \( \mathbf{C}^{n+1} \) does not change (as it should).
2. Assume that if the concentration at a specific time \( n \) in the sea is equal to \( \mathbf{C}_{eq} \) then the concentration stays constant at all later times
3. Add chemical reactions

In this chapter we will cover some theory related to the solution of nonlinear equations, and introduce the most used methods. A nonlinear problem is represented as a single equation or a system of equations, where the response is not changing proportionally to the input. Almost all physical systems are nonlinear, and one frequent use of the methods presented in this chapter is to determine model parameters by matching a nonlinear model to data.

Numerical methods that is guaranteed to find a solution (if it exists) are called closed methods, and open other vise. In many cases the closed methods requires more iterations for well behaved functions than the open methods. For one dimensional problems we will cover: fixed point iteration, bisection, Newton's method, and the secant method. For multidimensional problems we will cover Newton-Rapson method, which is a direct extension of Newton's method in one dimension, and the steepest decent. The main challenge is that there are (usually) more than one solution, the solution that you want for a specific problem is usually dictated by the underlying physics. If computational speed is not an issue, the method of choice is usually the bisection method. It is guaranteed to give an answer, but it might be slow. If speed is an issue, usually Newton's or the secant method will be the fastest (but it depends on the starting point). The secant method is sometimes preferred if the derivative of the function is costly to evaluate. Brents method is a method that combine the secant and bisection method (not covered), and is guaranteed to find a solution if the root is bracketed.

In many practical, engineering, applications one usually implements some of the methods described below directly inside functions. This is because it is usually faster than calling a separate all purpose nonlinear solver, and that one usually has a very good idea of what a good starting point for the nonlinear solver is.

Nonlinear equations

A nonlinear equation is simply an equation that is not linear. That means that when the variables changes the response is not changing proportional to the values of the variables. Solving a nonlinear equation always proceeds by iterations, we start with one or several initial guesses and then search for the solution. In many cases we do not know beforehand if the equation actually has a solution, or multiple solutions. An example of a nonlinear problem is:

Traditionally one collect all the terms on one side, to solve an equation of the form

In figure 22, the solution is shown graphically. Note that in one case the solution is when the graph of \( e^{-x} \), and \( x^2 \) intersect, whereas in the other case the root is located when \( x^2-e^{-x} \) intersect the $x-$axis.

In the case of more than one unknown, or a set of equations that must be satisfied simultaneously, equation \eqref{eq:nlin:fx} is replaced with a vector equation

Although this equation looks quite similar to equation \eqref{eq:nlin:fx}, this equation is much harder to solve. The only methods we will cover is the Newton Rapson method, which is a very good method if a good starting point is given. If you have a multidimensional problem, the advice is to try Newton-Raphson, if this method fails you need to try more advanced method, see e.g. [9].

Example: van der Waals equation of state

Before we begin with the numerical algorithms, let us consider an example: the van der Waals equation of state. The purpose is to illustrate some of the typical challenges. You are probably familiar with the ideal gas law:

where \( \nu=V/n \) is the molar volume of the gas, \( P \) is the pressure, \( V \) is the volume, \( T \) is the temperature, \( n \) is the number of moles of the gas, and \( R_g \) is the ideal gas constant. This equation is an example of an equation of state (EOS), it relates \( P \), \( T \), and \( \nu \). Thus if we know the pressure and temperature of the gas, we can calculate \( \nu \). Equation \eqref{eq:nlin:pvt} assumes that there are no interactions between the molecules in the gas. Clearly, this is too simplistic, and because of this one normally uses an EOS that better reflect the physical properties of the substance. A very famous EOS is the van der Waal EOS, which is a slight modification of equation \eqref{eq:nlin:pvt}:

\( a \) and \( b \) are material constants that needs to be determined experimentally. This equation is not used in industrial design, but most equations used in practice are based on equation \eqref{eq:nlin:vdw}. Multiplying equation \eqref{eq:nlin:vdw} with \( \nu^2 \), we get a non linear equation that is cubic in the molar volume. It turns out that cubic EOS are a class of equations that are quite successful in modeling the behavior of real systems [10]. However equation \eqref{eq:nlin:vdw} is a good starting point for more complex and realistic equations.

It is common practice to rescale EOS with respect to the critical point. At the critical point we have [ref]:

From equation \eqref{eq:nlin:crit1}, \eqref{eq:nlin:crit2}, and \eqref{eq:nlin:vdw}, it follows:

Inserting these equations into equation \eqref{eq:nlin:vdw}, and defining the reduced quantities \( \hat{P}=P/P_c \), \( \hat{T}=T/T_c \), \( \hat{\nu}=\nu/\nu_c \), we get

In figure 23, we have plotted the isotherms. Note that if \( \hat{T} < 1 \) (\( T < T_c \)), there might be more than one solution for the molar volume. This is clearly unphysical and additional constraints are needed. For the curve \( \hat{T}=0.9 \), the dashed lined shows that for \( \hat{P}=0.7 \), there are three solutions. This is a typical behavior of the cubic EOS, and physically it corresponds to the saturated case, where the vapor and liquid phase co-exist. The left root is the liquid state and the right root is the vapor state. The root in the middle represents a meta stable state.

Exercise 5.1: van der Waal EOS and CO$_2$

Use equation \eqref{eq:nlin:vdw}, and the parameters for CO$_2$: a=3.640 L$^2$bar/mol, and b=0.04267 L/mol, to test the van der Waal EOS in equation \eqref{eq:nlin:vdw}. Use that at 2 MPa and 100 $^\circ$C, CO$_2$ has a specific volume of 0.033586 m$^3$/kg.

Fixed-point iteration

A simple (but not always possible) way of solving a nonlinear equation is to reformulate the problem \( f(x)=0 \) to a problem of the form

The algorithm for solving this equation is to guess at a starting point, \( x_0 \), evaluate \( x_1=g(x_0) \), \( x_2=g(x_1) \), and so on. In some circumstances we might end up at a stable point, where \( x \) does not change. This point is termed a fixed point.

Note that the form of \( g(x) \) is not uniquely determined. For our function defined in equation \eqref{eq:nlin:exp}, we can solve for \( x \) directly

or we could write:

These functions are illustrated in figure 24, by visual inspection they look very similar, but as we will show in the next exercise the convergence is quite different.

Exercise 5.2: Implement the fixed point iteration

Write a Python function that utilizes the fixed point algorithm in the previous section, find the root of \( f(x)=x^2-e^{-x} \). In one case use \( g(x)=e^{-x/2} \), and in the other case use \( g(x)=x-x^2+e^{-x} \). How many iterations does it take in each case?

def iterative(x,g,prec=1e-8, MAXIT=1000):
    '''Approximate solution of x=g(x) by fixed point iterations.
    x : starting point for iterations 
    eps : desired precision
    Returns x when x does not change more than prec
    and number of iterations MAXIT are not exceeded
    '''
    eps = 1
    n=0
    while eps>prec and n < MAXIT:
        x_next = g(x)
        eps = np.abs(x-x_next)
        x = x_next
        n += 1
        if(np.isinf(x)):
            print('Quitting .. maybe bad starting point?')
            return x
    if (n<MAXIT):
        print('Found solution: ', x, ' After ', n, 'iterations')
    else:
        print('Max number of iterations exceeded')
    return x

Exercise 5.3: Finding the molar volume from the van der Waal EOS by fixed point iteration

Extend the code above to take as argument the van der Waal EOS. For simplicity we will use the rescaled EOS in equation \eqref{eq:nlin:vdwr}. Show that for the reduced temperature, $\hat{T}$=1.2, and pressure, $\hat{P}$=1.5, the reduced molar volume \( \hat{nu} \) is 1.3522091.

def dvdwEOS(nu,t,p):
    return (1+8*t/(p+3/nu**2))/3

def iterative(x,g,*args,prec=1e-8):
    MAX_ITER=1000
    eps = 1
    n=0
    while eps>prec and n < MAX_ITER:
        x_next = g(x,*args)
        eps = np.abs(x-x_next)
        x = x_next
        n += 1
    print('Number of iterations: ', n)
    return x

iterative(1,dvdwEOS,1.2,1.5)

When does the fixed point method fail?

If we replace \( e^{-x} \) with \( e^{1-x^2} \) in equation \eqref{eq:nlin:g3}, our method will not give a solution. You can easily verify that the \( x=1 \) is a solution, so why does our method fail? To investigate this in a bit more detail, we turn to Taylors formula (once again). Assume that the root is located at \( x^* \), and our guess is \( x_k \), then the next \( x \)-value will be

The true solution is \( x^* \), hence \( x^*=f(x^*) \), and we can write

where we have neglected higher order terms. The point is: at each iteration we want the distance \( x_1-x^* \) to decrease, i.e. to be smaller than \( x_0-x^* \). This can only be achieved if

In our example above we saw that if \( g(x)=x-x^2+e^{-x} \), we used 172 iterations and only 19 iterations if we replaced \( g(x) \) with \( h(x)=e^{-x/2} \) to converge to the same root $x$=0.70346742. We can now understand this, because \( g^\prime(x)=1-2x-e^{-x} \) and \( g(x^*)\simeq-0.90 \), whereas \( h^\prime(x)=-e^{-x/2}/2 \), and \( h^\prime(x^*)\simeq0.35 \). We expect the number of iterations, \( n \), needed to reach a certain precision, \( \varepsilon \), to scale as

We expect to use \( \log|h^\prime(x^*)|/\log|g^\prime(x^*)|\simeq10 \) more iterations using \( g(x) \) compared to \( h(x) \), which is close to the observed value of 172/19$\simeq 9$.

What to do when the fixed point method fails

taking the derivative with respect to x gives

Going from equation \eqref{eq:nlin:fpi1} to \eqref{eq:nlin:fpi2}, we have used the chain rule. Equation \eqref{eq:nlin:fpi3} is general, let us now specify to our fixed point iteration. Then we can use \( x^*=g(x^*)=y^* \), and \( x^*=g^{-1}(y^*)=g^{-1}(x^*) \) hence we can write the last equation as

Exercise 5.4: Solve \( x=e^{1-x^2} \) using fixed point iteration

The solution to \( x=e^{1-x^2} \) is clearly \( x=1 \).

• First try the fixed point method using \( g(x)=e^{1-x^2} \) to find the root \( x=1 \). Try to start very close to the true solution \( x=1 \). What is the value of \( g^\prime(x^*) \)?
• Next, invert \( g(x) \), what is the derivative of \( g^{-1}(x^*) \)? Try the fixed point method using \( g^{-1}(x^*) \)

Rate of convergence

The rate of convergence is the speed at which a convergent sequence approach the limit. Assume that our sequence \( x_{k} \) converges to the number \( x^* \), the sequence is said to converge linearly to \( x^* \) if there exists a number \( \mu\in < 0,1> \), such that

Inserting equation \eqref{eq:nlin:t2} in equation \eqref{eq:nlin:linconv}, we get:

Hence the fixed point iteration is expected to converge linearly to the correct solution. The definition in equation \eqref{eq:nlin:linconv}, can be extended to include the definition of quadratic, cubic, etc. convergence:

If \( q=2 \) the convergence is said to be quadratic and so on.

The bisection method

The idea behind bisection is that the root is bracketed, i.e. that there exists two points \( a \) and \( b \), such that \( f(a)\cdot f(b) < 0 \). In practice it might be a challenge to find these two points. However, if you know that the function has a only root between two values, and that speed is not a big issue this method guarantees that the root will be found within a finite number of steps. The basic idea behind the method is to divide the interval into two (i.e. bisecting the interval). The method only works if the function is continuous on the interval.

The algorithm is as follows:

• Test if \( f(a)\cdot f(b) < 0 \), if not return an error message
• Calculate the midpoint \( c=(a+b)/2 \). If \( f(a)\cdot f(c) < 0 \) the root is in the interval \( [a,c] \), else the root is in the interval \( [c,b] \)
• Half the interval, and test in which interval the root lies, and continue until a convergence criterion.

In figure 25, there is a graphical illustration. Below is an implementation of the bisection method.

def bisection(f,a,b,PREC=1e-8,MAXIT=100):
    '''Approximate solution of f(x)=0 on interval [a,b] by bisection.

    f   : f(x)=0.
    a,b : brackets the root f(a)*f(b) has to be negative 
    PREC: desired precision
    
    Returns the midpoint when it is closer than eps to the root, 
    unless MAXIT are not exceeded
    '''
    if f(a)*f(b) > 0:
        print('You need to bracket the root, f(a)*f(b) >= 0')
        return None
    prec=10*PREC
    c = 0.5*(a + b)
    for n in range(MAXIT):
        c_old = c 
        fc = f(c)
        if fc == 0:
            print('Found exact solution ', c, 
                    ' after ', n, 'iterations' )
            return c
        if f(a)*fc < 0:
            b = c
        else:
            a = c
        c = 0.5*(a+b)
        prec=np.abs(c_old-c)
    if n<MAXIT-1:
        print('Found solution ', c,' after ', n, 'iterations' )
        return c
    else:
        print('Max number of iterations: ', MAXIT, ' reached.') 
        print('Try to increase MAXIT')
        print('Returning best guess, value of function is: ', fc)
    return c

Rate of convergence

If \( c_n \) is the midpoint after \( n \) steps, the difference between the solution \( x^* \) and \( c_n \) is

Using our previous definition in equation \eqref{eq:nlin:qconv}, we find that

hence the bisection method converges linearly.

Newton's method

We can easily derive the algorithm by finding the formula for the tangent line. Using \( y=ax+b \) for the tangent line, we immediately know that \( a=f^\prime(x_k) \). \( b \) can be found as we know that the line intersects \( (x_k,f(x_k)) \): \( f(x_k)=f^\prime(x_k)x_k+b \), hence the equation for the tangent line is \( y=f^\prime(x_k)x+f(x_k)-f^\prime(x_k)x_k \). The next point is located where \( y \) crosses the \( x \)-axis, hence \( 0=f^\prime(x_k)x_{k+1}+f(x_k)-f^\prime(x_k)x_k \). Rearranging this equation, we can write Newtons method in the standard form

Note that the derivative of \( f(x) \) enters in equation \eqref{eq:nlin:newton}, which means that if our function has a extremal value in our search domain, Newtons method most likely will fail. In particular \( x_1 \), and \( x_4 \) in the figure to the right in figure 27 are bad starting point for Newtons method.

An implementation is shown below.

def newton(f,x, prec=1e-8,MAXIT=500):
    '''Approximate solution of f(x)=0 by Newtons method.
    The derivative of the function is calculated numerically
    f   : f(x)=0.
    x   : starting point  
    eps : desired precision
    
    Returns x when it is closer than eps to the root, 
    unless MAX_ITERATIONS are not exceeded
    '''
    MAX_ITERATIONS=MAXIT
    x_old = x
    h     = 1e-4
    for n in range(MAX_ITERATIONS):
        x_new = x_old - 2*h*f(x_old)/(f(x_old+h)-f(x_old-h))
        if(abs(x_new-x_old)<prec):
            print('Found solution:', x_new, 
                  ', after:', n, 'iterations.' )
            return x_new
        x_old=x_new
    print('Max number of iterations: ', MAXIT, ' reached.') 
    print('Try to increase MAXIT or decrease prec')
    print('Returning best guess, value of function is: ', f(x_new))
    return x_new

Comparing figure 25 and 26, you immediately get the sense that Newtons method converges faster, and indeed it does.

Rate of convergence

Newtons method is similar to the fixed point method, but where we do not use \( g(x)=x-f(x) \), but \( g(x)=x-\frac{f(x)}{f^\prime(x)} \). We will now analyze Newtons method, using the same approach as in the section When does the fixed point method fail?. First we expand \( g(x) \) around the root \( x^* \)

where we have skipped all higher order terms. You can easily verify that

\( x^* \) is a solution, hence \( f(x^*)=0 \), we then find from equation \eqref{eq:nlin:gn2} and \eqref{eq:nlin:gn3} that \( g^\prime(x^*)=0 \), and \( g^{\prime\prime}(x^*)=f^{\prime\prime}(x^*)/f^{\prime}(x^*)^2 \). Thus from equation \eqref{eq:nlin:nsec} we get

or equivalently:

The denominator has a power of two, and hence Newtons method is quadratic convergent (assuming that the sequence \( x_{k+1} \) is a convergent sequence). Note that it also follows from the analyses above that Newtons method will fail if the derivative at the root, \( f^\prime(x^*) \), is zero.

Exercise 5.5: Compare Newtons, Bisection and the Fixed Point method

Find the root of \( f(x)=x^2-e^{-x} \) using bisection, fixed point, and Newtons method, start at \( x=0 \). How many iterations do you need to use reach a precision of \( 10^{-8} \)? What happens if you widen the search domain or start further away from the root?

Secant method

The Newtons method is very good if you can choose a good starting point, and you can give in an analytical formula for the derivative. In some cases it is not possible to calculate the derivative analytically, then a very good method of choice is the secant method. It can be derived by simply replacing the derivative in Newtons method by the finite difference approximation

Inserting this equation into equation \eqref{eq:nlin:newton}, we get

For a graphical illustration see figure 29

Rate of convergence

The derivation of the rate of convergence for the secant method is a bit more involved. To simplify the notation we introduce the notation \( \varepsilon_k\equiv x_k-x^* \), where \( x^* \) is the exact solution. Subtracting \( x^* \) from each side of equation \eqref{eq:nlin:sec2} we get

we now make a Taylor expansion of \( f(x_k) \) and \( f(x_{k-1}) \) about the root \( x^* \)

where we have used the fact that \( f(x^*)=0 \). Inserting these equations into equation \eqref{eq:nlin:sec3} and neglecting terms of order \( \varepsilon_k^3 \) we get

where we have neglected higher powers of \( \varepsilon \). We are searching for a solution of the form \( \varepsilon_{k+1}=K\varepsilon_k^q \), \( q \) is the rate of convergence. We can invert this equation to get \( \varepsilon_k=K^{-1/q}\varepsilon_{k+1}^{1/q} \), or alternatively \( \varepsilon_{k-1}=K^{-1/q}\varepsilon_{k}^{1/q} \) (just set \( k\to k-1 \)). Inserting these equations into equation \eqref{eq:nlin:sec4}

Clearly, if this equation is to have a solution we must have

or \( q=1+1/q \). Solving this equation we get \( q=(1\pm\sqrt{5})/2 \), neglecting the negative solution, we find the rate of convergence for the secant method \( q=(1+\sqrt{5})/2\simeq 1.618 \).

Newton Rapson method

The derivation of Newtons method, equation \eqref{eq:nlin:newton}, done in the previous section was based on figure 26. We will now derive it using a slightly different approach, but which lends itself easier to extend Newtons method to higher dimensions. The starting point is to expand the function around \( x_k \), using Taylors formula

Equation \eqref{eq:nlin:newton} can be derived from equation \eqref{eq:nlin:nt} by simply demanding that we keep the linear terms, and that the next point \( x_{k+1} \) is located where the linear approximation intersects the \( x \)-axis, i.e. simply set \( f(x)=0 \), and \( x=x_{k+1} \) in equation \eqref{eq:nlin:nt}.

In higher order dimensions, we solve equation \eqref{eq:nlin:fvec}, and equation \eqref{eq:nlin:nt} is

\( \mathbf{J}(\mathbf{x}_k) \) is the Jacobian. As before, we simply set \( \mathbf{f}(\mathbf{x})=\mathbf{0} \), \( \mathbf{x}=\mathbf{x}_{k+1} \), and keep the linear terms, hence

To make the mathematics a bit more clear, let us specify to \( 2D \). Assume that \( \mathbf{f}(\mathbf{x})=[f_x(x,y),f_y(x,y)] \), then the Jacobian is

Gradient descent

This method used is to minimize functions (does not work for root finding). In many nonlinear problems, we would like to minimize (or maximize) a function. An ideal 2D example is shown in figure 30. The algorithm moves in the direction of steepest descent. Note that the step size might change towards the search.

Assume that we have a function \( \mathbf{f}(\mathbf{x}) \), that we would like to minimize. The gradient descent algorithm is simply to update parameters according to the derivative (gradient) of \( \mathbf{f} \)

\( \gamma \) is the learning rate, and a good choice of \( \gamma \) is important. \( \gamma \) might also change from one iteration to the other, and does not have to be constant.

Exercise 5.6: Gradient descent solution of linear regression

A very typical example is if we have a model and we would like to fit some parameters of the model to a data set (e.g. linear regression). Assume that we have observations \( (x_i,y_i) \) and model predictions \( f(x_i,\mathbf{\beta}) \), the model parameters are contained in the vector \( \mathbf{\beta} \). The least square, \( S \), is the square of the sum of all the residuals, i.e. the difference between the observations and model predictions

Specializing to linear regression, we choose the model to be linear

Equation \eqref{eq:nlin:lsq} now takes the form

The gradients are:

• Implement the gradient descent method using a constant learning rate of \( 10^{-3} \), to minimize the least square function
• Test the linear regression on the data set \( x_i=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9] \), and \( y=[1, 3, 2, 5, 7, 8, 8, 9, 10, 12] \), choose a starting value \( (b_0,b_1)=(0,0) \). What happens if you increase the learning rate?

def gradient_descent(f,x,df, g=.001, prec=1e-8,MAXIT=10):
    '''Minimize f(x) by gradient descent.
    f   : min(f(x))
    x   : starting point 
    df  : derivative of f(x)
    g   : learning rate
    prec: desired precision
    
    Returns x when it is closer than eps to the root, 
    unless MAXIT are not exceeded
    '''
    x_old = x
    for n in range(MAXIT):
        plot_regression_line(x_old)  
        x_new = x_old - g*df(x_old)
        if(abs(np.max(x_new-x_old))<prec):
            print('Found solution:', x_new, 
                  ', after:', n, 'iterations.' )
            return x_new
        x_old=x_new
    print('Max number of iterations: ', MAXIT, ' reached.') 
    print('Try to increase MAXIT or decrease prec')
    print('Returning best guess, value of function is: ', f(x_new))
    return x_new

x_obs_ = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) 
y_obs_ = np.array([1, 3, 2, 5, 7, 8, 8, 9, 10, 12]) 
def plot_regression_line(b,x=x_obs_, y=y_obs_): 
    global N_
    # plotting the actual points as scatter plot 
    plt.scatter(x, y, color = "m", 
               marker = "o", s = 30,label="data") 
  
    # predicted response vector 
    y_pred = b[0] + b[1]*x
  
    # plotting the regression line
    if(len(b)>1):
#        plt.plot(x, y_pred, color = "g", label = "R-squared = {0:.3f}".format(b[2]))
        plt.plot(x, y_pred, color = "g", label = "iteration:" + str(N_) +", (b[0],b[1])= ({0:.3f}".format(b[0]) + ", {0:.3f})".format(b[1]))
        plt.legend()
    else:
        plt.plot(x, y_pred, color = "g")
  
    # putting labels 
    plt.xlabel('x') 
    plt.ylabel('y') 
    plt.grid()
    plt.legend()
#    plt.savefig('../fig-nlin/stdec'+str(N_)+'.png', bbox_inches='tight',transparent=True)
    N_=N_+1  
    # function to show plot 
    plt.show() 


def Jacobian(x,f,dx=1e-5):
    N=len(x)
    x0=np.copy(x)
    f0=f(x)
    J=np.zeros(shape=(N,N))
    for j in range(N):
        x[j] = x[j] +  dx
        for i in range(N):   
            J[i][j] = (f(x)[i]-f0[i])/dx
        x[j] = x[j] -  dx
    return J




def newton_rapson(x,f,J=None, jacobian=False, prec=1e-8,MAXIT=100):
    '''Approximate solution of f(x)=0 by Newtons method.
    The derivative of the function is calculated numerically
    f   : f(x)=0.
    J   : Jacobian
    x   : starting point  
    eps : desired precision
    
    Returns x when it is closer than eps to the root, 
    unless MAX_ITERATIONS are not exceeded
    '''
    MAX_ITERATIONS=MAXIT
    x_old = np.copy(x)
    for n in range(MAX_ITERATIONS):
        plot_regression_line(x_old) 
        if not jacobian:
            J_=Jacobian(x_old,f)
        else:
            J_=J(x_old)
        z=np.linalg.solve(J_,-f(x_old))
        x_new=x_old+z
        if(np.sum(abs(x_new-x_old))<prec):
            print('Found solution:', x_new, 
                  ', after:', n, 'iterations.' )
            return x_new
        x_old=np.copy(x_new)
    print('Max number of iterations: ', MAXIT, ' reached.') 
    print('Try to increase MAXIT or decrease prec')
    print('Returning best guess, value of function is: ', f(x_new))
    return x_new


def gradient_descent(f,x,df, g=.001, prec=1e-8,MAXIT=10):
    '''Minimize f(x) by gradient descent.
    f   : min(f(x))
    x   : starting point 
    df  : derivative of f(x)
    g   : learning rate
    prec: desired precision
    
    Returns x when it is closer than eps to the root, 
    unless MAXIT are not exceeded
    '''
    x_old = x
    for n in range(MAXIT):
        plot_regression_line(x_old)  
        x_new = x_old - g*df(x_old)
        if(abs(np.max(x_new-x_old))<prec):
            print('Found solution:', x_new, 
                  ', after:', n, 'iterations.' )
            return x_new
        x_old=x_new
    print('Max number of iterations: ', MAXIT, ' reached.') 
    print('Try to increase MAXIT or decrease prec')
    print('Returning best guess, value of function is: ', f(x_new))
    return x_new
#end

def S(b,x=x_obs_,y=y_obs_):
    return np.sum((y-b[0]-b[1]*x)**2)

def dS(b,x=x_obs_,y=y_obs_):
    return np.array([-2*np.sum(y-b[0]-b[1]*x),
                     -2*np.sum((y-b[0]-b[1]*x)*x)])

def J(b,x=x_obs_,y=y_obs_):
    N=len(b)
    J=np.zeros(shape=(N,N))
    xs=np.sum(x)
    J[0][0]=2*len(x)
    J[0][1]=2*xs
    J[1][0]=2*xs
    J[1][1]=2*np.sum(x*x)
    return J
N_=0
print('Gradient ')
b=np.array([0,0])

Other useful methods

In this chapter we have covered the basic, but you should now be well equipped to dive into other methods. We highly recommend [9] as a starting point, although the code examples are written in C++, the theory is presented in a very accurate, but informal way.

• Brents method: uses root bracketing, bisection, and inverse quadratic interpolation. The 1D method of choice if the function and not its derivative is known

Algorithmic thinking

The only way to improve in coding and algorithmic thinking is practice. The concept of one dimensional numerical integration is easy to understand, i.e. to calculate the area under a curve. In this chapter we will implement several numerical methods, and it will serve as a very simple playground that illustrates the key aspects of numerical modeling

1. We start with a mathematical model (in this case an integral)
2. The mathematical model is formulated in discrete form
3. Then we design an algorithm to solve the model
4. The numerical solution for a test case is compared with the true solution (could be an analytical solution or data)
5. Error analysis: we investigate the accuracy of the algorithm by changing the number of iterations and/or make changes to the implementation or algorithm

The main point of this chapter is not to develop your own integration methods, the built in methods in Scipy will work in most cases. However, the way to break down the main task of calculating an integral into smaller tasks that is understandable by a computer, may work as a template for many different problems you would typically solve using a computer. A second motivation is that by analyzing the origin of numerical errors gives ideas for improving the algorithm, which is transferable to other problems.

The midpoint rule

Numerical integration is encountered in numerous applications in physics and engineering sciences. Let us first consider the most simple case, a function \( f(x) \), which is a function of one variable, \( x \). The most straight forward way of calculating the area \( \int_a^bf(x)dx \) is simply to divide the area under the function into \( N \) equal rectangular slices with size \( h=(b-a)/N \), as illustrated in figure 32. The area of one box is:

and the area of all the boxes is:

Note that the sum goes from \( k=0,1,\ldots,N-1 \), a total of \( N \) elements. We could have chosen to let the sum go from \( k=1,2,\ldots,N \). In Python, C, C++ and many other programming languages the arrays start by indexing the elements from \( 0,1,\ldots \) to \( N-1 \), therefore we choose the convention of having the first element to start at \( k=0 \).

Below is a Python code, where this algorithm is implemented for \( \int_0^\pi\sin (x)dx \)

import numpy as np
# Function to be integrated
def f(x):
    return np.sin(x)

def midpoint(f,a,b,N):
    """
    f : function to be integrated on the domain [a,b]
    N : number of integration points
    """
    h=(b-a)/N
    x=np.arange(a+0.5*h,b,h)
    return h*np.sum(f(x))
N=10
a=0
b=np.pi
Area = midpoint(f,a,b,N)
print('Numerical value= ', Area)
print('Error= ', (2-Area)) # Analytical result is 2

for k in range(0,N): # loop over k=0,1,..,N-1
    val = lower_limit+(k+0.5)*h # midpoint value
    area += func(val)
return area*h

The trapezoidal rule

The numerical error in the above example is quite low, only about 2$\%$ for \( N=5 \). However, by just looking at the graph above it seems likely that we can develop a better algorithm by using trapezoids instead of rectangles, see figure 33.

Earlier we approximated the area using the midpoint value: \( f(x_k+h/2)\cdot h \). Now we use \( A=A_1+A_2 \), where \( A_1=f(x_k)\cdot h \) and \( A_2=(f(x_k+h)-f(x_k))\cdot h/2 \), hence the area of one trapezoid is:

This is the trapezoidal rule, and for the whole interval we get:

Note that this formula was bit more involved to derive, but it requires only one more function evaluations compared to the midpoint rule. Below is a python implementation:

def trapezoidal(f,a,b,N):
    """
    f : function to be integrated on the domain [a,b]
    N : number of integration points
    """
    h=(b-a)/N
    x=np.arange(a+h,b,h)
    return h*(0.5*f(a)+0.5*f(b)+np.sum(f(x)))

In the table below, we have calculated the numerical error for various values of \( N \).

Note that we get the surprising result that this algorithm performs poorer, a factor of 2 than the midpoint rule. How can this be explained? By just looking at figure 32, we see that the midpoint rule actually over predicts the area from \( [x_k,x_k+h/2] \) and under predicts in the interval \( [x_k+h/2,x_{k+1}] \) or vice versa. The net effect is that for many cases the midpoint rule give a slightly better performance than the trapezoidal rule. In the next section we will investigate this more formally.

Numerical errors on integrals

It is important to know the accuracy of the methods we are using, otherwise we do not know if the computer produce correct results. In the previous examples we were able to estimate the error because we knew the analytical result. However, if we know the analytical result there is no reason to use the computer to calculate the result(!). Thus, we need a general method to estimate the error, and let the computer run until a desired accuracy is reached.

In order to analyze the midpoint rule in more detail we approximate the function by a Taylor series at the midpoint between \( x_k \) and \( x_k+h \):

Since \( f(x_k+h/2) \) and its derivatives are constants it is straight forward to integrate \( f(x) \):

The first term is simply the midpoint rule, to evaluate the two other terms we make the substitution: \( u=x-x_k \):

Note that all the odd terms cancels out, i.e \( \int_0^h(u-h/2)^m=0 \) for \( m=1,3,5\ldots \). Thus the error for the midpoint rule, \( E_{M,k} \), on this particular interval is:

where we have ignored higher order terms. We can easily sum up the error on all the intervals, but clearly \( f^{\prime\prime}(x_k+h/2) \) will not, in general, have the same value on all intervals. However, an upper bound for the error can be found by replacing \( f^{\prime\prime}(x_k+h/2) \) with the maximal value on the interval \( [a,b] \), \( f^{\prime\prime}(\eta) \):

where we have used \( h=(b-a)/N \). We can do the exact same analysis for the trapezoidal rule, but then we expand the function around \( x_k-h \) instead of the midpoint. The error term is then:

At the first glance it might look like the midpoint rule always is better than the trapezoidal rule, but note that the second derivative is evaluated in different points (\( \eta \) and \( \overline{\eta} \)). Thus it is possible to construct examples where the midpoint rule performs poorer than the trapezoidal rule.

Before we end this section we will rewrite the error terms in a more useful form as it is not so easy to evaluate \( f^{\prime\prime}(\eta) \) (since we do not know which value of \( \eta \) to use). By taking a closer look at equation \eqref{eq:numint:em}, we see that it is closely related to the midpoint rule for \( \int_a^bf^{\prime\prime}(x)dx \), hence:

The corresponding formula for the trapezoid formula is:

Practical estimation of errors on integrals (Richardson extrapolation)

From the example above we were able to estimate the number of steps needed to reach (at least) a certain precision. In many practical cases we do not deal with functions, but with data and it can be difficult to evaluate the derivative. We also saw from the example above that the algorithm gives a higher precision than what we asked for. How can we avoid doing too many iterations? A very simple solution to this question is to double the number of intervals until a desired accuracy is reached. The following analysis holds for both the trapezoid and midpoint method, because in both cases the (global) error scale as \( h^2 \) 1.

1: You can do the following analysis by assuming that the local error is \( h^3 \), but then you need to take into account that you need to take twice as many steps, which will give the same result.

Assume that we have evaluated the integral with a step size \( h_1 \), and the computed result is \( I_1 \). Then we know that the true integral is \( I=I_1+c h_1^2 \), where \( c \) is a constant that is unknown. If we now half the step size: \( h_2=h_1/2 \), then we get a new (better) estimate of the integral, \( I_2 \), which is related to the true integral \( I \) as: \( I=I_2+c h_2^2 \). Taking the difference between \( I_2 \) and \( I_1 \) give us an estimation of the error:

where we have used the fact that \( h_1=2h_2 \), Thus the error term is:

This might seem like we need to evaluate the integral twice as many times as needed. This is not the case, by choosing to exactly half the spacing we only need to evaluate for the values that lies halfway between the original points. We will demonstrate how to do this by using the trapezoidal rule, because it operates directly on the \( x_k \) values and not the midpoint values. The trapezoidal rule can now be written as:

in the last equation we have split the sum into odd an even values. The sum over the even values can be rewritten:

note that \( N_2 \) is replaced with \( N_1=N_2/2 \), we can now rewrite \( I_2 \) as:

Note that the first terms are actually the trapezoidal rule for \( I_1 \), hence:

The factor \( 1/2 \) in front of \( I_1(a,b) \), appears because \( h_2=h_1/2 \). A possible algorithm is then:

1. Choose a low number of steps to evaluate the integral, \( I_0 \), the first time, e.g. \( N_0=1 \)
2. Double the number of steps, \( N_1=2N_0 \)
3. Calculate the missing values by summing over the odd number of steps \( \sum_{k=\text{odd values}}^{N_1-1}f(a+k h_1) \)
4. Check if \( E_1(a,b)=\frac{1}{3}(I_1-I_0) \) is lower than a specific tolerance
5. If yes quit, if not, return to 2, and continue until \( E_i(a,b)=\frac{1}{3}(I_{i+1}-I_{i}) \) is lower than the tolerance

Below is a Python implementation:

def int_adaptive_trapez2(lower_limit, upper_limit,func,tol):
    """
    adaptive quadrature, integrate a function from lower_limit
    to upper_limit within tol*(upper_limit-lower_limit)

    """
    S=[]
    S.append([lower_limit,upper_limit])
    I=0
    iterations=0
    while S:
        iterations +=1
        a,b=S.pop(-1) # last element
        m=(b+a)*0.5   # midpoint
        I1=0.5*(b-a)*(func(a)+func(b)) #trapezoidal for 1 interval 
        I2=0.25*(b-a)*(func(a)+func(b)+2*func(m)) #trapezoidal for 2 intervals
        if(np.abs(I1-I2)<3*np.abs((b-a)*tol)):
            I+=I2     # accuarcy met
        else:
            S.append([a,m]) # half the interval 
            S.append([m,b])
    print("Number of iterations: ", iterations)
    return I

If you compare the number of terms used in the adaptive trapezoidal rule, which was developed by halving the step size, and the adaptive midpoint rule that was derived on the basis of the theoretical error term, you will find the adaptive midpoint rule is more efficient. So why go through all this trouble? In the next section we will see that the development we did for the adaptive trapezoidal rule is closely related to Romberg integration, which is much more effective.

Romberg integration

The adaptive algorithm for the trapezoidal rule in the previous section can be easily improved by remembering that the true integral was given by 2 : \( I=I_i+ch_i^2+\mathcal{O}(h^4) \). The error term was in the previous example only used to check if the desired tolerance was achieved, but we could also have added it to our estimate of the integral to reach an accuracy to fourth order:

2: Note that all odd powers of \( h \) is equal to zero, thus the corrections are always in even powers.

As before the error term \( \mathcal{O}(h^4) \), can be written as: \( ch^4 \). Now we can proceed as in the previous section: First we estimate the integral by one step size \( I_i=I+ch_i^4 \), next we half the step size \( I_{i+1}=I+ch_{i+1}^4 \) and use these two estimates to calculate the error term:

but now we are in the exact situation as before, we have not only the error term but the correction up to order \( h^4 \) for this integral:

Each time we half the step size we also gain a higher order accuracy in our numerical algorithm. Thus, there are two iterations going on at the same time; one is the iteration that half the step size (\( i \)), and the other one is the increasing number of higher order terms added (which we will denote \( m \)). We need to improve our notation, and replace the approximation of the integral (\( I_i \)) with \( R_{i,m} \). Equation \eqref{eq:numint:rom}, can now be written:

A general formula valid for any \( m \) can be found by realizing:

where, as before \( h_{i-1}=2h_i \). Subtracting equation \eqref{eq:numint:rom0} and \eqref{eq:numint:rom1}, we find an expression for the error term:

Then the estimate for the integral in equation \eqref{eq:numint:rom1} is:

A possible algorithm is then:

1. Evaluate \( R_{0,0}=\frac{1}{2}\left[f(a)+f(b)\right](b-a) \) as the first estimate
2. Double the number of steps, \( N_{i+1}=2N_i \) or half the step size \( h_{i+1}=h_i/2 \)
3. Calculate the missing values by summing over the odd number of steps \( \sum_{k=\text{odd values}}^{N_1-1}f(a+k h_{i+1}) \)
4. Correct the estimate by adding all the higher order error term \( R_{i,m+1}=R_{i,m}+\frac{1}{4^m-1}(R_{i+1,m+1}-R_{i,m+1}) \)
5. Check if the error term is lower than a specific tolerance \( E_{i,m}(a,b)=\frac{1}{4^{m+1}-1}(R_{i,m}-R_{i-1,m}) \), if yes quit, if no goto 2, increase \( i \) and \( m \) by one

The algorithm is illustrated in figure 34.

Note that the tolerance term is not the correct one as it uses the error estimate for the current step, which we also use correct the integral in the current step to reach a higher accuracy. Thus the error on the integral will always be lower than the user specified tolerance. Below is a Python implementation:

def int_romberg(func,a, b,tol,show=False):
    """ calculates the area of func on the domain [a,b]
        for the given tol, if show=True the triangular
        array of intermediate results are printed """
    Nmax = 100
    R = np.empty([Nmax,Nmax]) # storage buffer
    h = (b-a) # step size
    R[0,0]    =.5*(func(a)+func(b))*h
    N = 1
    for i in range(1,Nmax):
        h /= 2
        N *= 2
        odd_terms=0
        for k in range (1,N,2): # 1, 3, 5, ... , N-1
            val        = a + k*h
            odd_terms += func(val)
		# add the odd terms to the previous estimate	
        R[i,0]   = 0.5*R[i-1,0] + h*odd_terms 
        for m in range(i): 
			# add all higher order terms in h
            R[i,m+1]   = R[i,m] + (R[i,m]-R[i-1,m])/(4**(m+1)-1)                  
		# check tolerance, best guess			
        calc_tol = abs(R[i,i]-R[i-1,i-1])       
        if(calc_tol<tol):
            break  # estimated precision reached 
    if(i == Nmax-1):
        print('Romberg routine did not converge after ',
              Nmax, 'iterations!')
    else:      
        print('Number of intervals = ', N)

    if(show==True):
        elem = [2**idx for idx in range(i+1)]
        print("Steps StepSize Results")
        for idx in range(i+1):
            print(elem[idx],' ',
                  "{:.6f}".format((b-a)/2**idx),end = ' ')
            for l in range(idx+1):
                print("{:.6f}".format(R[idx,l]),end = ' ')
            print('')  
    return R[i,i] #return the best estimate

Note that the Romberg integration only uses 32 function evaluations to reach a precision of \( 10^{-8} \), whereas the adaptive midpoint and trapezoidal rule in the previous section uses 20480 and 9069 function evaluations, respectively.

Alternative implementation of adaptive integration

Before we proceed, we will consider an alternative implementation of the adaptive method presented in the previous sections, with the following modification

1. We will use Simpsons rule (see the exercise at the end), which takes the following form \( \int_a^bf(x)dx\simeq\frac{h}{6}\left[f(a)+4f(a+\frac{h}{2})+2f(a+h)+ 4f(a+3\frac{h}{2})+2f(a+2h)+\cdots+f(b)\right] \)
2. We only divide the intervals needed to reach the desired accuracy.

Simpsons rule is accurate up to \( \mathcal{O}(h^4) \), and by following the same arguments as above we can estimate the error as \( E_i(a,b)=\frac{1}{15}(I_{i+1}-I_{i}) \). The factor 1/15 (as opposed to 1/3) originates from the higher order accuracy. The integration proceeds as follows

• S is an empty list
• S.append([a,b])
• \( I=0 \)
• while S not empty do:
• [a,b]=S.pop(-1)
• \( m=(b+a)/2 \)
• \( I_1= \) simpson_step(a,b)
• \( I_2= \) simpson_step(a,m)+simpson_step(m,b)
• if \( |I_1-I_2| < 15|b-a|\cdot tol \)
• \( I+=I_2 \)
• else:
• S.append([a,m])
• S.append([m,b])
• return \( I \)

Note the use of the list S, we remove the interval \( [a,b] \) from the list and calculates the integral. If the integral is not accurate enough we add to new intervals to the list, and continue until we reach the desired accuracy, then we proceed with the next interval. Since we remove (pop) the element from the list, we know that we will finish the evaluation once the list is empty. This algorithm allows for different sub interval to have different degrees of subdivisions, contrary to Rombergs algorithm. The full python implementation is shown below

def simpson_step(a, b,func):
    m=0.5*(a+b)
    return (b-a)/6*(func(a)+func(b)+4*func(m))

def int_adaptive_simpson(func,a, b,tol):
    """
    adaptive quadrature, integrate a function from a
    to b within tol*(b-a) uses simpsons rule
    """
    S=[]
    S.append([a,b])
    I=0
    iterations=0
    while S:
        iterations +=1
        a,b=S.pop(-1) # last element
        m=(b+a)*0.5   # midpoint
        I1=simpson_step(a,b,func) #simpsons for 1 interval 
        I2=simpson_step(a,m,func)+simpson_step(m,b,func) # ...2 intervals
        if(np.abs(I1-I2)<15*np.abs((b-a)*tol)):
            I+=I2     # accuarcy met
        else:
            S.append([a,m]) # half the interval 
            S.append([m,b])
    print("Number of iterations: ", iterations)
    return I

Gaussian quadrature

Many of the methods we have looked into are of the type:

where the function is evaluated at fixed interval. For the midpoint rule \( \omega_k=h \) for all values of \( k \), for the trapezoid rule \( \omega_k=h/2 \) for the endpoints and \( h \) for all the interior points. For the Simpsons rule (see exercise) \( \omega_k=h/3, 4h/3,2h/3,4h/3,\ldots,4h/3,h/3 \). Note that all the methods we have looked at so far samples the function in equal spaced points, \( f(a+k h) \), for \( k=0, 1, 2\ldots, N-1 \). If we now allow for the function to be evaluated at unevenly spaced points, we can do a lot better. This realization is the basis for Gaussian Quadrature. We will explore this in the following, but to make the development easier and less cumbersome, we transform the integral from the domain \( [a,b] \) to \( [-1,1] \):

The factor in front comes from the fact that \( dt=(b-a)dx/2 \), thus we can develop our algorithms on the domain \( [-1,1] \), and then do the transformation back using: \( t=(b-a)x/2+(b+a)/2 \).

This idea is best understood by a couple of examples. Assume that we want to use \( N=1 \) in equation \eqref{eq:numint:qq1}:

We now choose \( f(x) \) to be a polynomial of as large a degree as possible, but with the requirement that the integral is exact. If \( f(x)=1 \), we get:

hence \( \omega_0=2 \). If we choose \( f(x)=x \), we get:

hence \( x_0=0 \).

This equation is equal to the midpoint rule, by choosing \( b=a+h \) we reproduce equation \eqref{eq:numint:mid0}. If we choose \( N=2 \):

we can show that now $ f(x)=1,\,x,\,x^2\,x^3$ can be integrated exact:

hence there are four unknowns and four equations. The solution is: \( \omega_0=\omega_1=1 \) and \( x_0=-x_1=1/\sqrt{3} \).

def int_gaussquad2(func, lower_limit, upper_limit):
    x   = np.array([-1/np.sqrt(3.),1/np.sqrt(3)])
    w   = np.array([1, 1])
    xp  = 0.5*(upper_limit-lower_limit)*x
    xp += 0.5*(upper_limit+lower_limit)
    area = np.sum(w*func(xp))
    return area*0.5*(upper_limit-lower_limit)

The case N=3

For the case \( N=3 \), we find that \( f(x)=1,x,x^2,x^3,x^4,x^5 \) can be integrated exactly:

the solution to these equations are \( \omega_{0,1,2}=5/9, 8/9, 5/9 \) and \( x_{1,2,3}=-\sqrt{3/5},0,\sqrt{3/5} \). Below is a Python implementation:

def int_gaussquad3(lower_limit, upper_limit,func):
    x  = np.array([-np.sqrt(3./5.),0.,np.sqrt(3./5.)])
    w  = np.array([5./9., 8./9., 5./9.])
    xp = 0.5*(upper_limit-lower_limit)*x
    xp += 0.5*(upper_limit+lower_limit)
    area = np.sum(w*func(xp))
    return area*0.5*(upper_limit-lower_limit)

Note that the Gaussian quadrature converges very fast. From \( N=2 \) to \( N=3 \) function evaluation we reduce the error (in this specific case) from 6.5% to 0.1%. Our standard trapezoidal formula needs more than 20 function evaluations to achieve this, the Romberg method uses 4-5 function evaluations. How can this be? If we use the standard Taylor formula for the function to be integrated, we know that for \( N=2 \) the Taylor formula must be integrated up to \( x^3 \), so the error term is proportional to \( h^4f^{(4)}(\xi) \) (where \( \xi \) is some x-value in \( [a,b] \)). \( h \) is the step size, and we can replace it with \( h\sim (b-a)/N \), thus the error scale as \( c_N/N^4 \) (where \( c_N \) is a constant). Following the same argument, we find for \( N=3 \) that the error term is \( h^6f^{(6)}(\xi) \) or that the error term scale as \( c_N/N^6 \). Each time we increase \( N \) by a factor of one, the error term reduces by \( N^2 \). Thus if we evaluate the integral for \( N=10 \), increasing to \( N=11 \) will reduce the error by a factor of \( 11^2=121 \).

Error term on Gaussian integration

The Gaussian integration rule of order \( N \) integrates exactly a polynomial of order \( 2N-1 \). From Taylors error formula, see equation \eqref{eq:taylor:error} in the chapter Finite differences, we can easily see that the error term must be of order \( 2N \), and be proportional to \( f^{(2N)}(\eta) \), see [7] for more details on the derivation of error terms. The drawback with an analytical error term derived from series expansion is that it involves the derivative of the function. As we have already explained, this is very unpractical and it is much more practical to use the methods described in the section sec:numint:parct. Let us consider this in more detail, assume that we evaluate the integral using first a Gaussian integration rule with \( N \) points, and then \( N+1 \) points. Our estimates of the "exact" integral, \( I \), would then be:

In principle \( h_{N+1}\neq h_{N} \), but in the following we will assume that \( h_N\simeq h_{N+1} \), and \( h\ll 1 \). Subtracting equation \eqref{eq:numint:gerr1} and \eqref{eq:numint:gerr2} we can show that a reasonable estimate for the error term \( ch^{2N} \) would be:

If this estimate is lower than a given tolerance we can be quite confident that the higher order estimate \( I_{N+1} \) approximate the true integral within our error estimate. This is the method implemented in SciPy, integrate.quadrature

Common weight functions for classical Gaussian quadratures

Integrating functions over an infinite range

Integrating a function over an infinite range can be done by the following trick. Assume that we would like to evaluate

If we introduce the following substitution

then if \( x=a \), \( z=0 \), and if \( x\to\infty \) then \( z\to1 \), hence:

Exercise 6.1: Numerical Integration

a) Show that for a linear function, \( y=a\cdot x+b \) both the trapezoidal rule and the rectangular rule are exact

b) Consider \( I(a,b)=\int_a^bf(x)dx \) for \( f(x)=x^2 \). The analytical result is \( I(a,b)=\frac{b^3-a^3}{3} \). Use the Trapezoidal and Midpoint rule to evaluate these integrals and show that the error for the Trapezoidal rule is exactly twice as big as the Midpoint rule.

c) Use the fact that the error term on the trapezoidal rule is twice as big as the midpoint rule to derive Simpsons formula: \( I(a,b)=\sum_{k=0}^{N-1}I(x_k,x_k+h)=\frac{h}{6}\left[f(a)+ 4f(a+\frac{h}{2})+2f(a+h)+4f(a+3\frac{h}{2})+2f(a+2h)+\cdots+f(b)\right] \) Hint: \( I(x_k,x_k+h)=M(x_k,x_k+h)+E_M \) (midpoint rule) and \( I(x_k,x_k+h)=T(x_k,x_k+h)+E_T=T(x_k,x_k+h)-2E_M \) (trapezoidal rule).

d) Show that for \( N=2 \) (\( f(x)=1,x,x^3 \)), the points and Gaussian quadrature rule for \( \int_{0}^{1}x^{1/2}f(x)\,dx \) is \( \omega_{0,1}=-\sqrt{70}{150} + 1/3, \sqrt{70}{150} + 1/3 \) and \( x_{0,1}=-2\sqrt{70}{63} + 5/9, 2\sqrt{70}{63} + 5/9 \)

1. Integrate \( \int_0^1x^{1/2}\cos x\,dx \) using the rule derived in the exercise above and compare with the standard Gaussian quadrature rule for (\( N=2 \), and \( N=3 \)).

e) Make a Python program that uses the Midpoint rule to integrate experimental data that are unevenly spaced and given in the form of two arrays.

0D models

By 0D models we mean models that are only dependent on one variable, usually time. These models are often called "lumped parameter models", because they ignore spatial dependency and the parameters in the model are adjusted to match experimental data. In many cases it is more important to predict how a system evolves in time, than the spatial dependency. A very recent example is the spread of infectious deceases, like Covid-19, where the evolution of total number of infected people might be the most important and not where the infection happens.

Ordinary differential equations

Physical systems evolves in space and time, and very often they are described by a ordinary differential equations (ODE) and/or partial differential equations (PDE). The difference between an ODE and a PDE is that an ODE only describes the changes in one spatial dimension or time, whereas a PDE describes a system that evolves in the \( x- \), \( y- \), \( z- \) dimension and/or in time. In the following we will spend a significant amount of time to explore one of the simplest algorithm, Eulers method. Sometimes this is exactly the algorithm you would like to use, but with very little extra effort much more sophisticated algorithms can easily be implemented, such as the Runge-Kutta fourth order method. However, all these algorithms, will at some point run into the same kind of troubles if used reckless. Thus we will use the Eulers method as a playground, investigate when the algorithm run into trouble and suggests ways to fix it, these approaches can easily be extended to the higher order methods. Most of the other algorithms boils down to the same idea of extrapolating a function using derivatives multiplied with a small step size.

A simple model for fluid flow

Let us consider a simple example from chemical engineering, a continuous stirred tank reactor (CSTR), see figure 35. The flow is incompressible (\( q_\text{out}=q_\text{in} \)), a fluid is entering on the top and exiting at the bottom, the tank has a fixed volume \( V \). Assume that the tank is filled with saltwater, and that freshwater is pumped into it, how much time does it take before \( 90\% \) of the saltwater is replaced with freshwater? The tank is well mixed, illustrated with the propeller, this means that at every time the concentration is uniform in the tank, i.e. that \( C(t)=C_\text{out}(t) \).

The concentration \( C \) is measured in gram of salt per liter water, and the flow rate \( q \) is liter of water per day. The model for the salt balance in this system can be described in words by:

In our case there are no generation of salt within the system so this term is zero. The flow of salt into the system during a time \( \Delta t \) is: \( q_\text{in}(t)\cdot C_\text{in}(t)\cdot \Delta t=q(t)\cdot C_\text{in}(t)\cdot \Delta t \), the flow of salt out of the system is: \( q_\text{out}(t)\cdot C_\text{out}(t)\cdot \Delta t=q(t)\cdot C(t)\cdot \Delta t \), and the accumulation during a time step is: \( C(t+\Delta t)\cdot V - C(t)\cdot V \), hence:

Note that it is not a priori apparent, which time the concentrations and flow rates on the right hand side should be evaluated at, we could have chosen to evaluate them at \( t+\Delta t \), or at any time \( t\in [t,t+\Delta t] \). We will return to this point later in this chapter. Dividing by \( \Delta t \), and taking the limit \( \Delta t\to 0 \), we can write equation \eqref{eq:ode:cstr1} as:

Seawater contains about 35 gram salt/liter fluid, if we assume that the fresh water contains no salt, we have the boundary conditions \( C_\text{in}(t)=0 \), $C(0)=$35gram/l. The equation \eqref{eq:ode:cstr2} the reduces to:

this equation can easily be solved, by dividing by \( C \), multiplying by \( dt \) and integrating:

This equation can be inverted to give \( t=-\tau\ln[C(t)/C] \). If we assume that the volume of the tank is 1m$^3$=1000liters, and that the flow rate is 1 liter/min, we find that $\tau$=1000min=0.69days and that it takes about $-0.69\ln0.9\simeq1.6$days to reduce the concentration by 90$\%$ to 3.5 gram/liter.

Euler's method

If the system gets slightly more complicated, e.g several tanks in series with a varying flow rate or if salt was generated in the tank, there is a good chance that we have to solve the equations numerically to obtain a solution. Actually, we have already developed a numerical algorithm to solve equation \eqref{eq:ode:cstr2}, before we arrived at equation \eqref{eq:ode:cstr2} in equation \eqref{eq:ode:cstr1}. This is a special case of Eulers method, which is basically to replace the derivative in equation \eqref{eq:ode:cstr2}, with \( (C(t+\Delta t)-C(t))/\Delta t \). By rewriting equation \eqref{eq:ode:cstr1}, so that we keep everything related to the new time step, \( t+\Delta t \), on one side, we get:

we introduce the short hand notation: \( C(t)=C_n \), and \( C(t+\Delta t)=C_{n+1} \), hence the algorithm can be written more compact as:

In the script below, we have implemented equation \eqref{eq:ode:eu2}.

def analytical(x):
    return np.exp(-x)

def euler_step(c_old, c_in, tau_inv,dt):
    fact=dt*tau_inv
    return (1-fact)*c_old+fact*c_in

def ode_solv(c_into,c_init,t_final,vol,q,dt):
    f=[];t=[]
    tau_inv = q/vol
    c_in    = c_into #freshwater into tank
    c_old   = c_init #seawater present 
    ti=0.
    while(ti <= t_final):
        t.append(ti); f.append(c_old)
        c_new = euler_step(c_old,c_in,tau_inv,dt)     
        c_old = c_new
        ti   += dt
    return t,f

In figure 36 the result of the implementation is shown for different values of \( \Delta t \). Clearly we see that the results are dependent on the step size, as the step increases the numerical solution deviates from the analytical solution. At some point the numerical algorithm fails completely, and produces results that have no meaning.

Error analysis - Euler's method

There are two obvious questions:

1. When does the algorithm produce unphysical results?
2. What is an appropriate step size?

Let us consider the first question, clearly when the concentrations gets negative the solution is unphysical. From equation \eqref{eq:ode:eu2}, we see that when \( \Delta t/\tau > 1 \), the concentration become negative. For this specific case (the CSTR), there is a clear physical interpretation of this condition. Inserting \( \tau=V/q \), we can rewrite the condition \( \Delta t/\tau < 1 \) as \( q\Delta t < V \). The volume into the tank during one time step is: \( q\Delta t \), which means that whenever we flush more than one tank volume through the tank during one time step, the algorithm fails. When this happens the new concentration in the tank cannot be predicted from the old one. This makes sense, because we could have switched to a new solution (e.g. seawater) during that time step, then the new solution does not have any relation to the old solution.

The second question, "what is an appropriate step size?", is a bit more difficult to answer. One strategy could be to simply use the results from chapter [Taylor], where we showed that the truncation error had a minimum value with a step size of \( 10^{-8} \) (when using a first order Taylor approximation). How does the value \( 10^{-8} \) relate to the step sizes in minutes used in our Euler implementation? In order to see the connection, we need to rewrite equation \eqref{eq:ode:cstr2} in a dimensionless form, by making the following substitution: \( t\to t/\tau \):

As we found earlier $\tau = 1000$min, thus a step size of e.g. 1 min would correspond to a dimensionless time step of $\Delta t\to$1min/1000min$=10^{-3}$. This number can be directly compared to the value \( 10^{-8} \), which is the lowest value we can choose without getting into trouble with round off errors on the machine.

As already mentioned a step size of \( 10^{-8} \), is probably the smallest we can choose with respect to round off errors, but it is smaller than necessary and would lead to large simulation times. If it takes 1 second to run the simulation with a step size of \( 10^{-3} \), it would take \( 10^5 \) seconds or 1 day with a step size of \( 10^{-8} \). To continue the error analyses, we write our ODE for a general system as:

or in discrete form:

\( h \) is now the (dimensionless) step size, equal to \( \Delta t \) if the derivative is with respect to \( t \) or \( \Delta x \) if the derivative is respect to \( x \) etc. Note that we have also included the error term related to the numerical derivative, \( \eta_n\in[t_n,t_n+h] \). At each step we get an error term, and the distance between the true solution and our estimate, the local error, after \( N \) steps is:

Note that when we replace the sum with an integral in the equation above, this is only correct if the step size is not too large. From equation \eqref{eq:ode:eu3} we see that even if the error term on the numerical derivative is \( h^2 \), the local error is proportional to \( h \) (one order lower). This is because we accumulate errors for each step.

In the following we specialize to the CSTR, to see if we can gain some additional insight. First we change variables in equation \eqref{eq:ode:cstr3}: \( y=C(t)/C_0 \), and \( x=t/\tau \), hence:

The solution to this equation is \( y(x)=e^{-x} \), substituting back for the new variables \( y \) and \( x \), we reproduce the result in equation \eqref{eq:ode:sol}. The local error, equation \eqref{eq:ode:eu3}, reduces to:

we have assumed that \( x_0=t_0/\tau=0 \). This gives the estimated local error at time \( x_f \). For \( x_f=0 \), the numerical error is zero, this makes sense because at \( x=0 \) we know the exact solution because of the initial conditions. When we move further away from the initial conditions, the numerical error increases, but equation \eqref{eq:ode:eu4} ensures us that as long as the step size is low enough we can get as close as possible to the true solution, since the error scales as \( h \) (at some point we might run into trouble with round off error in the computer).

Can we prove directly that we get the analytical result? In this case it is fairly simple, if we use Eulers method on equation \eqref{eq:ode:simple}, we get:

or alternatively:

In the last equation, we have used the the fact the number of steps, \( N \), is equal to the simulation time divided by the step size, hence: \( N=x_f/h \). From calculus, the equation above is one of the well known limits for the exponential function: \( \lim_{x\to\infty}(1+k/x)^{mx}=e^{mk} \), hence:

when \( h\to0 \). Below is an implementation of the Euler algorithm in this simple case, we also estimate the local error, and global error after \( N \) steps.

import matplotlib.pyplot as plt
import numpy as np
def euler(tf,h):
    t=[];f=[]
    ti=0.;fi=1.
    t.append(ti);f.append(fi)
    global_err=0.
    while(ti<= tf):
        ti+=h
        fi=fi*(1-h)
        global_err += abs(np.exp(-ti)-fi)
        t.append(ti);f.append(fi)
    print("error= ", np.exp(-ti)-fi," est.err=", .5*h*(1-np.exp(-ti)))
    print("global error=",global_err)
    return t,f
                                        
t,f=euler(1,1e-5)

By changing the step size \( h \), you can easily verify that the local error systematically increases or decreases proportional to \( h \). Something curious happens with the global error when the step size is changed, it does not change very much. The global error involves a second sum over the local error for each step, which can be approximated as a second integration in equation \eqref{eq:ode:eu4}:

Note that the global error does not go to zero when the step size decreases, which can easily be verified by changing the step size. This is strange, but can be understood by the following argument: when the step size decreases the local error scales as \( \sim h \), but the number of steps scales as \( 1/h \), so the global error must scale as \( h\times 1/h \) or some constant value. Usually it is much easier to control the local error than the global error, this should be kept in mind if you ever encounter a problem where it is important control the global error. For the higher order methods that we will discuss later in this chapter, the global error will go to zero when \( h \) decreases.

The answer to our original question, ''What is an appropriate step size?'', will depend on what you want to achieve in terms of local or global error. In most practical situations you would specify a local error that is acceptable for the problem under investigation and then choose a step size where the local error always is lower than this value. In the next subsection we will investigate how to achieve this in practice.

Adaptive step size - Euler's method

We want to be sure that we use a step size that achieves a certain accuracy in our numerical solution, but at the same time that we do not waste simulation time using a too low step size. The following approach is similar to the one we derived for the Romberg integration, and a special case of what is known as Richardson Extrapolation. The method is easily extended to higher order methods.

We know that Eulers algorithm is accurate to second order. Our estimate of the new value, \( y_1^* \) (where we have used a$\,{}^*$ to indicate that we have used a step size of size \( h \)), should then be related to the true solution \( y(t_1) \) in the following way:

The constant \( c \) is unknown, but it can be found by taking two smaller steps of size \( h/2 \). If the steps are not too large, our new estimate of the value \( y_1 \) will be related to the true solution as:

The factor 2 in front of \( c \) is because we now need to take two steps, and we accumulate a total error of \( 2c(h/2)^2=ch^2/2 \). It might not be completely obvious that the constant \( c \) should be the same in equation \eqref{eq:ode:aeb0} and \eqref{eq:ode:aeb1}. If you are not convinced, there is an exercise at the end of the chapter. We define:

The truncation error in equation \eqref{eq:ode:aeb1} is:

Now we have everything we need: We want the local error to be smaller than some predefined tolerance, \( \epsilon^\prime \), or equivalently that \( \epsilon\le\epsilon^\prime \). To achieve this we need to use an optimal step size, \( h^\prime \), that gives us exactly the desired error:

Dividing equation \eqref{eq:ode:ae6} by equation \eqref{eq:ode:ae5b}, we can estimate the optimal step size:

where the estimated error, \( \epsilon \), is calculated from equation \eqref{eq:ode:ae5b}. Equation \eqref{eq:ode:ae7} serves two purposes, if the estimated error \( \epsilon \) is higher than the tolerance, \( \epsilon^\prime \), we have specified it will give us an estimate for the step size we should choose in order to achieve a higher accuracy, if on the other hand \( \epsilon^\prime > \epsilon \), then we get an estimate for the next, larger step. Before the implementation we note, as we did for the Romberg integration, that equation \eqref{eq:ode:ae5b} also gives us an estimate for the error term in equation \eqref{eq:ode:aeb1} as an improved estimate of \( y_1 \). This we get for free and will make our Euler algorithm accurate to \( h^3 \), hence the improved Euler step, \( \hat{y_1} \), is to subtract the error term from our previous estimate:

Below is an implementation of the adaptive Euler algorithm:

def one_step(c_old, c_in,h):
    return (1-h)*c_old+h*c_in

def adaptive_euler(c_into,c_init,t_final,tol=1e-4):
    f=[];t=[]
    c_in    = c_into #freshwater into tank
    c_old   = c_init #seawater present 
    ti=0.; h_new=1e-3;
    no_steps=0
    global_err=0.
    while(ti <= t_final):
        t.append(ti); f.append(c_old)
        toli=10.*tol; # a high init tolerance to enter while loop
        while(toli>tol):# first two small steps
            hi=h_new
            k1 = one_step(c_old,c_in,hi*.5)
            k2 = one_step(k1,c_in,hi*.5)
            # ... and one large step
            k3 = one_step(c_old,c_in,hi)
            toli = abs(k3-k2)
            h_new=hi*np.sqrt(tol/toli)
            no_steps+=3
        toli=1.
        c_old=2*k2-k3 # higher order correction
 # normal Euler, uncomment and inspect the global error
 #       c_old = k2 
        ti   += hi
        global_err += abs(np.exp(-ti)-c_old)
    print("No steps=", no_steps, "Global Error=", global_err)
    return t,f

In figure 37 the result of the implementation is shown. Note that the number of steps for an accuracy of \( 10^{-6} \) is only about 3000. Without knowing anything about the accuracy, we would have to assume that we needed a step size of the order of \( h \) in order to reach a local accuracy of \( h \) because of equation \eqref{eq:ode:eu3}. In the current case, we would have needed \( 10^7 \) steps, which would lead to unnecessary long simulation times.

Runge-Kutta methods

The Euler method only have an accuracy of order \( h \), and a global error that do not go to zero as the step size decrease. The Runge-Kutta methods may be motivated by inspecting the Euler method in figure 38. The Euler method uses information from the previous time step to estimate the value at the new time step. The Runge Kutta methods uses the information about the slope between the points \( t_n \) and \( t_n+h \). By inspecting figure 38, we clearly see that by using the slope at \( t_n+h/2 \) would give us a significant improvement. The 2. order Runge-Kutta method can be derived by Taylor expanding the solution around \( t_n+h/2 \), we do this by setting \( t_n+h=t_n+h/2+h/2 \):

Similarly we can expand the solution in \( y(t_n) \) about \( t_n+h/2 \), by setting \( t_n=t_n+h/2-h/2 \):

Subtracting these two equations the term \( y(t_n+\frac{h}{2}) \), and all even powers in the derivative cancels out:

In the last equation, we have used equation \eqref{eq:ode:ode}. Note that we now have an expression that is very similar to Eulers algorithm, but it is accurate to order \( h^3 \). There is one problem, and that is that the function \( f \) is to be evaluated at the point \( y_{n+1/2}=y(t_n+h/2) \) which we do not know. This can be fixed by using Eulers algorithm: \( y_{n+1/2}=y_n+h/2f(y_n,t_n) \). We can do this even if Eulers algorithm has an error term of order \( h^2 \), because the \( f \) in equation \eqref{eq:ode:rk3} is multiplied by \( h \), and thus our algorithm is still has an error term of order \( h^3 \).

Below is a Python implementation of equation \eqref{eq:ode:rk4}:

def fm(c_old,c_in):
    return c_in-c_old

def rk2_step(c_old, c_in, h):
    k1=h*fm(c_old,c_in)
    k2=h*fm(c_old+0.5*k1,c_in)
    return c_old+k2

def ode_solv(c_into,c_init,t_final,h):
    f=[];t=[]
    c_in  = c_into #freshwater into tank
    c_old = c_init #seawater present 
    ti=0.
    while(ti <= t_final):
        t.append(ti); f.append(c_old)
        c_new = rk2_step(c_old,c_in,h)     
        c_old = c_new
        ti   += h
    return t,f

In figure 39 the result of the implementation is shown. Note that when comparing Runge-Kutta 2. order with Eulers method, see figure 39 and 36, we of course have the obvious result that a larger step size can be taken, without loosing numerical accuracy. It is also worth noting that we can take steps that is larger than the tank volume. Eulers method failed whenever the time step was larger than one tank volume (\( h=t/\tau>1 \)), whereas the Runge-Kutta method finds a physical solution for step sizes lower than twice the tank volume. If the step size is larger, we see that the concentration in the tank increases, which is clearly unphysical.

The Runge-Kutta fourth order method is one of he most used methods, it is accurate to order \( h^4 \), and has an error of order \( h^5 \). The development of the algorithm itself is similar to the 2. order method, but of course more involved. We just quote the result:

In figure 40 the result of the Runge-Kutta fourth order is shown, by comparing it to figure 39 it is easy to see that a larger step size can be chosen.

Adaptive step size - Runge-Kutta method

Just as we did with Eulers method, we can implement an adaptive method. The derivation is exactly the same, but this time our method is accurate to fourth order, hence the error term is of order \( h^5 \). We start by taking one large step of size \( h \), our estimate, \( y_1^* \) is related to the true solution, \( y(t_1) \), in the following way:

Next, we take two steps of half the size, \( h/2 \), hence:

Subtracting equation \eqref{eq:ode:rka0} and \eqref{eq:ode:rka1}, we find an expression similar to equation \eqref{eq:ode:ae5}:

or \( c=16\Delta/(15h^5) \). For the Euler scheme, \( \Delta \) also happened to be equal to the truncation error, but in this case it is:

we want the local error, \( \epsilon \), to be smaller than some tolerance, \( \epsilon^\prime \). The optimal step size, \( h^\prime \), that gives us exactly the desired error is then:

Dividing equation \eqref{eq:ode:rka3} by equation \eqref{eq:ode:rka5}, we can estimate the optimal step size:

\( \epsilon \) can be calculated from equation \eqref{eq:ode:rka5}. In figure 41 the result of an implementation is shown (see the exercises).

In general we can use the same procedure any method accurate to order \( h^p \), and you can easily verify that:

Conservation of mass

A mathematical model of a physical system should always be formulated in such a way that it is consistent with the laws of nature. In practical situations this statement is usually equivalent to state that the mathematical model should respect conservation laws. The conservation laws can be conservation of mass, energy, momentum, electrical charge, etc. In our example with the mixing tank, we were able to derive an expression for the concentration of salt out of the tank, equation \eqref{eq:ode:sol}, by demanding conservation of mass (see equation \eqref{eq:ode:cstr1}).

A natural question to ask is then: If our mathematical model respect conservation of mass, are we sure that our solution method respect conservation of mass? We of course expect that when the grid spacing approaches zero our numerical solution will get closer and closer to the analytical solution. Clearly when \( \Delta x\to 0 \), the mass is conserved. So what is the problem? The problem is that in many practical problems we cannot always have a step size that is small enough to ensure that our solution always is close enough to the analytical solution. The physical system we consider might be very complicated (e.g. a model for the earth climate), and our ODE system could be a very small part of a very big system. A very good test of any code is to investigate if the code respect the conservation laws. If we know that our implementation respect e.g. mass conservation at the discrete level, we can easily test mass conservation by summing up all the mass entering, and subtracting the mass out of and present in our system. If the mass is not conserved exactly, there is a good chance that there is a bug in our implementation.

If we now turn to our system, we know that the total amount of salt in the system when we start is \( C(0)V \). The amount entering is zero, and the amount leaving each time step is \( q(t)C(t)\Delta t \). Thus we should expect that if we add the amount of salt in the tank to the amount that has left the system we should always get an amount that is equal to the original amount. Alternatively, we expect \( \int_{t_0}^t qC(t)dt + C(t)V -C(0)V=0 \). Adding the following code in the while(ti <= t_final): loop:

mout += 0.5*(c_old+c_new)*q*dt
mbal  = (c_new*vol+mout-vol*c_init)/(vol*c_init)

it is possible to calculate the amount of mass lost (note that we have used the trapezoidal formula to calculate the integral). In the table below the fraction of mass lost relative to the original amount is shown for the various numerical methods.

We clearly see from the table that the Runge-Kutta methods performs better than Eulers method, but all of the methods violates mass balance.

This might not be a surprise as we know that our numerical solution is always an approximation to the analytical solution. How can we then formulate an algorithm that will respect conservation laws at the discrete level? It turns out that for Eulers method it is not so difficult. Eulers algorithm at the discrete level (see equation \eqref{eq:ode:eu0}) is actually a two-step process: first we inject the fresh water while we remove the ``old`` fluid and then we mix. By thinking about the problem this way, it makes more sense to calculate the mass out of the tank as \( \sum_kq_kC_k\Delta t_k \). If we in our implementation calculates the mass out of the tank as:

mout += c_old*q*dt
mbal  = (c_new*vol+mout-vol*c_init)/(vol*c_init)

We easily find that the mass is exactly conserved at every time for Eulers method. The concentration in the tank will of course not be any closer to the analytical solution, but if our mixing tank was part of a much bigger system we could make sure that the mass would always be conserved if we make sure that the mass out of the tank and into the next part of the system was equal to \( qC(t)\Delta t \).

Solving a set of ODE equations

What happens if we have more than one equation that needs to be solved? If we continue with our current example, we might be interested in what would happen if we had multiple tanks in series. This could be a very simple model to describe the cleaning of a salty lake by injecting fresh water into it, but at the same time this lake was connected to two nearby fresh water lakes, as illustrated in figure 42. The weakest part of the model is the assumption about complete mixing, in a practical situation we could enforce complete mixing with the salty water in the first tank by injecting fresh water at multiple point in the lake. For the two next lakes, the degree of mixing is not obvious, but salt water is heavier than fresh water and therefore it would sink and mix with the fresh water. Thus if the flow rate was slow, one might imaging that a more or less complete mixing could occur. Our model then could answer questions like, how long time would it take before most of the salt water is removed from the first lake, and how much time would it take before most of the salt water was cleared from the whole system? The answer to these questions would give practical input on how much and how fast one should inject the fresh water to clean up the system. If we had data from an actual system, we could compare our model predictions with data from the physical system, and investigate if our model description was correct.

For simplicity we will assume that all the lakes have the same volume, \( V \). The governing equations follows as before, by assuming mass balance (equation \eqref{eq:ode:mbal}):

Taking the limit \( \Delta t\to 0 \), we can write equation \eqref{eq:ode:cstr3a} as:

Let us first derive the analytical solution: Only the first tank is filled with salt water \( C_0(0)=C_{0,0} \), \( C_1(0)=C_2(0)=0 \), and \( C_\text{in}=0 \). The solution to equation \eqref{eq:ode:cstr3b} is, as before \( C_0(t)=C_{0,0}e^{-t/\tau} \), inserting this equation into equation \eqref{eq:ode:cstr3c} we find:

where we have use the technique of integrating factors when going from equation \eqref{eq:ode:cstr3e} to \eqref{eq:ode:cstr3f}. Inserting equation \eqref{eq:ode:cstr3g} into equation \eqref{eq:ode:cstr3d}, solving the equation in a similar way as for \( C_1 \) we find:

The numerical solution follows the exact same pattern as before if we introduce a vector notation. Before doing that, we rescale the time \( t\to t/\tau \) and the concentrations, \( \hat{C_i}=C_i/C_{0,0} \) for \( i=0,1,2 \), hence:

In figure 43 results of an implementation using Runge-Kutta 4. order is shown (see exercises for more details).

Stiff sets of ODE and implicit methods

As already mentioned a couple of times, our system could be part of a much larger system. To illustrate this, let us now assume that we have two tanks in series. The first tank is similar to our original tank, but the second tank is a sampling tank, 1000 times smaller.

The governing equations can be found by requiring mass balance for each of the tanks (see equation \eqref{eq:ode:mbal}:

Taking the limit \( \Delta t\to 0 \), we can write equation \eqref{eq:ode:cstr2a} as:

Assume that the first tank is filled with seawater, \( C_0(0)=C_{0,0} \), and fresh water is flooded into the tank, i.e. \( C_\text{in}=0 \). Before we start to consider a numerical solution, let us first find the analytical solution: As before the solution for the first tank (equation \eqref{eq:ode:cstr2bb}) is:

where \( \tau_0\equiv V_0/q \). Inserting this equation into equation \eqref{eq:ode:cstr2b}, we get:

where \( \tau_1\equiv V_1/q \).

Next, we will consider the numerical solution. You might think that these equations are more simple to solve numerically than the equations with three tanks in series discussed in the previous section. Actually, this system is much harder to solve with the methods we have discussed so far. The reason is that there are now two time scales in the system, \( \tau_1 \) and \( \tau_2 \). The smaller tank sets a strong limitation on the step size we can use, because we should never use step sizes larger than a tank volume. Thus if you use the code in the previous section to solve equation \eqref{eq:ode:cstr2bb} and \eqref{eq:ode:cstr2b}, it will not find the correct solution, unless the step size is lower than \( 10^{-3} \). Equations of this type are known as stiff.

These types of equations are often encountered in practical applications. If our sampling tank was extremely small, maybe \( 10^6 \) smaller than the chemical reactor, then we would need a step size of the order of \( 10^{-8} \) or lower to solve the system. This step size is so low that we easily run into trouble with round off errors in the computer. In addition the simulation time is extremely long. How do we deal with this problem? The solution is actually quite simple. The reason we run into trouble is that we require that the concentration leaving the tank must be a small perturbation of the old one. This is not necessary, and it is best illustrated with Eulers method. As explained earlier Eulers method can be viewed as a two step process: first we inject a volume (and remove an equal amount: \( qC(t)\Delta t \)), and then we mix. Clearly when we try to remove more than what is left, we run into trouble. What we want to do is to remove or flood much more than one tank volume through the tank during one time step, this can be achieved by \( q(t)C(t)\Delta t\to q(t+\Delta t)C(t+\Delta t)\Delta t \). The term \( q(t+\Delta t)C(t+\Delta t)\Delta t \) now represents the mass out of the system during the time step \( \Delta t \).

The methods we have considered so far are known as explicit, whenever we replace the solution in the right hand side of our algorithm with \( y(t+\Delta t) \) or (\( y_{n+1} \)), the method is known as implicit. Implicit methods are always stable, meaning that we can take as large a time step that we would like, without getting oscillating solution. It does not mean that we will get a more accurate solution, actually explicit methods are usually more accurate.

Let us consider our example further, and for simplicity use the implicit Eulers method:

This equation is equal to equation \eqref{eq:ode:cstr2a}, but the concentrations on the right hand side are now evaluated at the next time step. The immediate problem is now that we have to find an expression for \( C_{n+1} \) that is given in terms of known variables. In most cases one needs to use a root finding method, like Newtons method, in order to solve equation \eqref{eq:ode:cstr2ai}. In this case it is straight forward to show:

In figure 45 the result of the implementation is shown, note that quite large step sizes can be used without inducing non physical results.

Exercise 7.1: Truncation error in Euler's method

In the following we will take a closer look at the adaptive Eulers algorithm and show that the constant \( c \) is indeed the same in equation \eqref{eq:ode:aeb0} and \eqref{eq:ode:aeb1}. The true solution \( y(t) \), obeys the following equation:

and Eulers method to get from \( y_0 \) to \( y_1 \) by taking one (large) step, \( h \) is:

We will also assume (for simplicity) that in our starting point \( t=t_0 \), the numerical solution, \( y_0 \), is equal to the true solution, \( y(t_0) \), hence \( y(t_0)=y_0 \).

a) Show that when we take one step of size \( h \) from \( t_0 \) to \( t_1=t_0+h \), \( c=y^{\prime\prime}(t_0)/2 \) in equation \eqref{eq:ode:aeb0}.

b) Show that when we take two steps of size \( h/2 \) from \( t_0 \) to \( t_1=t_0+h \), Eulers algorithm is:

c) Find an expression for the local error when using two steps of size \( h/2 \), and show that the local error is: \( \frac{1}{2}ch^2 \)

Monte Carlo methods

Monte Carlo methods are named after the Monte Carlo Casino in Monaco, this is because at its core it uses random numbers to solve problems. Monte Carlo methods are quite easy to program, and they are usually much more intuitive than a theoretical approach.

Monte Carlo integration ''hit and miss''

Let us start with a simple illustration of the Monte Carlo Method (MCM), Monte Carlo integration. To the left in figure 46 there is a shape of a pond. Imagine that we wanted to estimate the area of the pond, how could we do it? Assume further that you did not have your phone or any other electronic devices to help you.

One possible approach is: First to walk around it, and put up some bands (illustrated by the black dotted line). Then estimate the area inside the bands (e.g. 4$\times$3 meters). Then we would know that the area was less than e.g. 12m$^2$. Finally, and this is the difficult part, throw rocks randomly inside the bands. The number of rocks hitting the pond divided by the total number rocks thrown should be equal to the area of the pond divided by the total area inside the bands, i.e. the area of the pond should be equal to:

It is important that we throw the rocks randomly, otherwise equation \eqref{eq:mc:mci} is not correct. Now, let us investigate this in more detail, and use the idea of throwing rocks to estimate \( \pi \). To the right in figure 46, there is a well known shape, a circle. The area of the circle is \( \pi d^2/4 \), and the shape is given by \( x^2+y^2=d^2/4 \). Assume that the circle is inscribed in a square with sides of \( d \). To throw rocks randomly inside the square, is equivalent pick random numbers with coordinates \( (x,y) \), where \( x\in[0,d] \) and \( y\in[0,d] \). We want all the \( x- \) and $y-$values to be chosen with equal probability, which is equivalent to pick random numbers from a uniform distribution. Below is a Python implementation:

import numpy as np

def estimate_pi(N,d):
#   random.seed(2)
    D2=d*d/4; dc=0.5*d
    A=0
    for k in range(0,N):
        x=np.random.uniform(0,d)
        y=np.random.uniform(0,d)
        if((x-dc)**2+(y-dc)**2 <= D2):
            A+=1
    # estimate area of circle: d*d*A/N
    return 4*A/N

N=1000;d=2
pi_est=estimate_pi(N,d)
print('Estimate for pi= ', pi_est,' Error=', np.pi-pi_est)

In the table below, we have run the code for \( d=1 \) and different values of \( N \).

We clearly see that a fair amount of rocks or numbers needs to be used in order to get a good estimate. If you run this code several times you will see that the results changes from time to time. This makes sense as the coordinates \( x \) and \( y \) are chosen at random.

for k in range(0,N):
        x=np.random.uniform(0,d)
        y=np.random.uniform(0,d)

x=np.random.uniform(0,d,size=N)
y=np.random.uniform(0,d,size=N)

Random number generators

There are much to be said about random number generators. The MCM depends on a good random number generator, otherwise we cannot use the results from statistics to develop our algorithms. Below, we briefly summarize some important points that you should be aware of:

1. Random number generators are generally of two types: hardware random number generator (HRNG) or pseudo random number generator (PRNG).
2. HRNG uses a physical process to generate random numbers, this could atmospheric noise, radioactive decay, microscopic fluctuations, which is translated to an electrical signal. The electrical signal is converted to a digital number (1 or 0), by sampling the random signal random numbers can be generated. The HRNG are often named true random number generators, and their main use are in cryptography.
3. PRNG uses a mathematical algorithm to generate an (apparent) random sequence. The algorithm uses an initial number, or a seed, to start the sequence of random number. The sequence is deterministic, and it will generate the same sequence of numbers if the same seed is used. At some point the algorithm will reproduce itself, i.e. it will have certain period. For some seeds the period may be much shorter.
4. Many of the PRNG are not considered to be cryptographically secure, because if a sufficiently long sequence of random numbers are generated from them, the rest of the sequence can be predicted.
5. Python uses the Mersenne Twister algorithm to generate random numbers, and has a period of \( 2^{19937}−1\simeq4.3\cdot10^{6001} \). It is not considered to be cryptographically secure.

In Pythons random.uniform function, a random seed is chosen each time the code is run, but if we set e.g. random.seed(2), the code will generate the same sequence of numbers each time it is called.

Encryption

This section can be skipped as it is not relevant for development of the numerical algorithms, but it is a good place to explain the basic idea behind encryption of messages. A very simple, but not a very good encryption, is to replace all the letters in the alphabet with a number, e.g. A=1, B=2, C=3, etc. This is what is know as a substitution cipher, it does not need to be a number it could be a letter, a sequence of letters, letters and numbers etc. The receiver can solve the code by doing the reverse operation.

The weakness of this approach is that it can fairly easily be cracked, by the following approach: First we analyze the encrypted message and find the frequency of each of the symbols. Assume that we know that the message is written in English, then the frequency of symbols can be compared with the frequency of letters from a known English text (the most common is E (12$\%$), then T (9$\%$), etc.). We would then guess that the most occurring symbol probably is an E or T. When some of the letters are in place, we can compare with the frequency of words, and so on. By the help of computers this process can easily be automated.

A much better algorithm is to not replace a letter with the same symbol. To make it more clear, consider our simple example where A=1, B=2, C=3, \( \ldots \). If we know say that A=1 but we add a random number, then our code would be much harder to crack. Then the letter A could be several places in the message but represented as a complete different number. Thus we could not use the frequency of the various symbols to crack the message.

How can the receiver decrypt the message? Obviously, it can be done if both the sender and receiver have the same sequence of random numbers (or the key). This can be achieved quite simple with random number generators, if we know the seed used we can generate the same sequence of random numbers. If Alice where to send a message to Bob without Eve knowing what it is, Alice and Bob could agree to send a message that was scrambled using Pythons Mersenne-Twister algorithm with seed=2.

The weakness of this approach is of course that Eve could convince Alice or Bob to give her the seed or the key. Another possibility is that Eve could write a program that tested different random number generators and seeds to decipher the message. How to avoid this?

Let us assume that Alice and Bob each had their own hardware random generator. This generator generated random numbers that was truly random, and the sequence could not be guessed by any outsider. Alice do not want to share her key (sequence of random numbers) with Bob, and Bob would not share his key with Alice. How can they send a message without sharing the key? One possible way of doing it is as follows: Alice write a message and encrypt it with her key, she send the message to Bob. Bob then encrypt the message with his key, he sends it back to Alice. Alice then decrypt the message with her key and send it back to Bob. Now, Bob can decrypt it with his own key and read the message. The whole process can be visualized by thinking of the message as box with the message. Alice but her padlock on the box (keeps her key for her self), she sends the message to Bob. Bob locks the box with his padlock, now there are two padlocks on the box. He sends the box back to Alice, Alice unlocks her padlock with her key, and sends it back to Bob. The box now only has Bob's key, he can unlock the box and read the message. The important point is that the box was never unlocked throughout the transaction, and Alice and Bob never had to share the key with anyone.

Errors on Monte Carlo integration and the binomial distribution

How many rocks do we need to throw in order to reach a certain accuracy? To answer this question we need some results from statistics. Our problem of calculating the integral is closely related to the binomial distribution. When we throw a rock one of two things can happen i) the rock falls into the water, or ii) it falls outside the pond. If we denote the probability that the rock falls into the pond as \( p \), then the probability that it falls outside the pond, \( q \), has to be \( q=1-p \). This is simply because there are no other possibilities and the sum of the two probabilities has to be one: \( p+q=p+(1-p)=1 \). The binomial distribution is given by:

\( p(k) \) is the probability that an event happens \( k \) times after \( n \) trials. The mean, \( \mu \), and the variance, \( \sigma^2 \), of the binomial distribution is:

Before we proceed, we should take a moment and look a little more into the meaning of equation \eqref{eq:mc:bin} to appreciate its usefulness. A classical example of the use of the binomial formula is to toss a coin, if the coin is fair it will have an equal probability of giving us a head or tail, hence \( p=1/2 \). Equation \eqref{eq:mc:bin}, can answer questions like: ''What is the probability to get only heads after 4 tosses?''. Let us calculate this answer using equation \eqref{eq:mc:bin}, the number of tosses is 4, the number of success is 4 (only heads each time)

''What is the probability to get three heads in four tosses?'', using the same equation, we find:

In figure 47, all the possibilities are shown. The number of possibilities are 16, and there are only one possibility that we get only heads, i.e. the probability is 1/16 as calculated in equation \eqref{eq:mc:coin}. In the figure we also see that there are 4 possible ways we can get three heads, hence the probability is 4/16=1/4 as calculated in equation \eqref{eq:mc:coin2}.

Now, let us return to our original question, ''What is the error on our estimate of the integral, when using the MCM?''. Before we continue we should also clean up our notation, let \( I \) be the value of the true integral, \( A \) is our estimate of the integral, and \( I_N \) is the area of the rectangle. First, let us show that the mean or expectation value of the binomial distribution is related to our estimate of the area of the pond or the circle, \( A \). In our case we draw \( n=N \) random numbers, and \( k \) times the coordinate falls inside the circle, equation \eqref{eq:mc:binm} tells us that the mean value is \( np \). \( p \) is the probability that the coordinate is within the area to be integrated, hence as before \( p \) is equal to the area to be integrated divided by the area of the total domain, thus:

or

Equation \eqref{eq:mc:bins}, gives us an estimate of the variance of the mean value. Assume for simplicity that we can replace \( 1-p\simeq p \), this is of course only correct if the area of the rectangle is twice as big as our pond, but we are only interested in an estimate of the error, hence \( \sigma^2\simeq np^2 \). We can now use the standard deviation as an estimate of the error of our integral:

In the last equation we have replaced \( p \) with \( A/I_N \). Hence, the error of our integral is inversely proportional to the square root of the number of points.

The mean value method

How does our previous method compare with some of our standard methods, like the midpoint rule? The error for the MC method scales as \( 1/\sqrt{N} \), in our previous error estimates we used the step length, \( h \), as an indicator of the accuracy, and not \( N \). The s$N$ is related to the number of points as \( h=(b-a)/n \), where \( b \) and \( a \) are the integration limit. Thus our MCM scales as \( 1/\sqrt{n}\sim h^{1/2} \), this is actually worse than the midpoint or trapezoidal rule, which scaled as \( h \).

The MCM can be improved. We will first describe the mean value method. In the last section we calculated the area of a circle by picking random numbers inside a square and estimated the fraction of points inside the circle. This is equivalent to calculate the area of a half circle, and multiply with 2:

The half-circle is now centered at the origin. Before we proceed we write our integral in a general form as:

Instead of counting the number of points inside the curve given by \( f(x) \), we could instead use the mean of the function, which we will define as \( \overline{f}=\sum_k f(x_k)/N \):

Note that this formula is similar to the midpoint rule, but now the function is not evaluated at the midpoint, but at several points and we use the average value.

Below is an implementation:

import numpy as np

def f(x,D2):
    return 2*np.sqrt((D2-x*x))

def mcm_mean(N,d):
    D2=d*d/4
    x=np.random.uniform(-d/2,d/2,size=N)
    A=np.sum(f(x,D2))
    # estimate for area: A/N
    return d*A/N

In the table below we have compared the mean value method with the ''hit and miss'' method. We see that the mean value method performs somewhat better, but there are some random fluctuations and in some cases it performs poorer.

We also see that in this case the error scales as \( 1/\sqrt{N} \).

Basic properties of probability distributions

The MCM is closely tied to statistics, and it is important to have a basic understanding of probability density functions (PDF). In the previous section, we used a random number generator to give us random numbers in an interval. All the numbers are picked with an equal probability. Another way to state this is to say that: we draw random numbers from an uniform distribution. Thus all the numbers are drawn with an equal probability \( p \). What is the value of \( p \)? That value is given from another property of PDF's, all PDF's must be normalized to 1. This is equivalent to state that the sum of all probabilities must be equal to one. Thus for a general PDF, \( p(x) \), we must have:

A uniform distribution, \( p(x)=U(x) \), is given by:

you can easily verify that \( \int_{-\infty}^{\infty}U(x)=1 \). In the MCM we typically evaluate expectation values. The expectation value, \( E[f] \), for a function is defined:

specializing to a uniform distribution, \( p(x)=U(x) \), we get:

Rearranging this equation, we see that we can write the above equation as:

This equation is the same as equation \eqref{eq:mc:I2}, but in the previous section we never explained why the expectation value of \( f(x_k) \) was equal to the integral. The derivation above shows that \( \int_a^bf(x)dx \) is equal to the expectation value of \( f(x) \) only under the condition that we pick numbers from a uniform distribution.

To make this a bit more clearer, let us specialize to \( f(x)=x \). In this case the expectation value is equal to the mean:

In this case we also find that the numerical error scales as \( N^{-1/2} \), from the definition of the variance

Example: Monte Carlo integration of a hyper sphere

The volume of a hyper sphere is known:

where \( D \) is the number of dimensions \( \Gamma(D/2+1) \) is the gamma function, if \( n \) is an integer then \( \Gamma(n)=(n-1)! \) and \( \Gamma(n+1/2)=(2n)!/(4^nn!)\sqrt{\pi} \). You can easily verify that for \( D=2,3 \), \( V(R)=\pi R^2, 4/3\pi R^3 \), respectively.

In the case of MC integration, we simply place the sphere inside a cube, and then count the number of points that hits inside the hyper sphere:

def mc_nballII(N=1000,D=3,R=1):
    """
    Calculates the volume of a hypersphere using accept and reject
    N: Number of random points
    D: Number of dimensions
    R: Radius of hypersphere
    """
    r=0
    for d in range(D):
        xi=np.random.uniform(-R,R,size=N)
        r+=xi*xi
    r=np.sqrt(r)
    vol = np.sum(r<=R)
    print("Number of points inside n-ball: ", vol)
    return vol/N*(2*R)**D