The aim of the Lab is to learn to use finite difference approximations of derivatives of a function and to derive finite difference methods for boundary value problems of ordinary differential equations.
Let us start from the numerical differentiation by finite difference approximations. Often we have a situation, where we can compute values of certain function (like option pricing function) but can not compute exact derivatives, since we do not have analytic formulas for the function (or the formulas are so complicated that computing the formula for a derivative is too much work). Then we can use finite difference formulas which allow us to compute approximate values of derivatives by using only values of the function of interest at certain points.
Well-known finite difference approximations are as follows: $$f'(a)\approx \frac{f(a+h)-f(a-h)}{2h}\qquad \text{(error}\sim ch^2),$$ $$f''(a)\approx \frac{f(a-h)-2f(a)+f(a+h)}{h^2}\qquad \text{(error}\sim ch^2).$$ These formulas tell us, how we should compute the approximate value and also how the error goes to 0 when $h$ goes to 0. The fact that the error behaves like $c\,h^2$ for an unknown constant $c$ when $h$ tends to 0 means that we usually get very accurate results even for relatively large values of $h$. We usually do not know the value of the constant $c$ (it depends on the function we are considering), but knowing the behavior of the error is the key for some practical ways of estimating the actual error of obtained results.
Write a function myDerivative that takes four arguments - a name of a function, the value of $a$, the value of $h$ and the order of the derivative (1 or 2) - and computes the value of the specified derivative at $x$ using the finite difference approximation given above. Using this function, verify the accuracy of the error estimates in the case of $f(x)=e^x$, $f(x)=\sqrt{x}$ at $a=1$: compute the derivatives for $h=1,\frac{1}{2},\frac{1}{4},\ldots,\frac{1}{2^{10}}$ and find the quotient of the error to $h^2$ for both the first and second order derivatives. The quotients should approach a constant value.
Solution
import numpy as np
from scipy import linalg
##Exercise 1
def myDerivative(f,a,h,order):
if order==1:
answer=???
return answer
Let us check that the function works correctly:
a=0
print(myDerivative(np.exp,a,0.1,1)-np.exp(a))#should be 0.0016675001984409743
print(myDerivative(np.exp,a,0.1,2)-np.exp(a)) #should be 0.0008336111607227803
Now check if the error of the approximate derivatives computed by myDerivative behaves like $c\cdot h^2$ for some constant c for small enough $h$:
#define a vector containing h values
h=??
#compute and print the quotients (errors divided by h**2)
As we can see, the quotients tend to a constant value in all cases.
Theoretically the smaller value of $h$ we use, the better approximation to the derivative we get and previous exercise showed that this is indeed the case for $h$ values down to $2^{-10}=\frac{1}{1024}$. But theoretical error estimates do not take into account that computers use only a fixed number of bytes to store a real number and therefore all numbers can not be stored exactly. Actually, all numbers are rounded so that approximately 15 correct digits are stored. Therefore, if we compute the difference of two numbers which are very close to each other, only a few correct digits of the result (except leading zeros) may be found. For example, if actual numbers were $A=1.21111111$ and $B=1.20000001$ but we could store only 3 digits, the stored numbers would be $a=1.21$ and $b=1.20$. If we use $a$ and $b$ for computing $A-B$, we get only 1 digit (after leading zeros) right - we get 0.01, but correct difference is 0.0111111. This means that the relative error of the answer is quite large and if we multiply the result with a large number (like 1000, for example), the error of the result may be quite significant. This may happen when we use finite difference formulas for computing derivatives in the case of small $h$, since $f(a+h)$ and $f(a-h)$ may be quite close to each other and their difference may have a large relative error. Let us see that this is indeed the case.
Compute the errors of the finite difference approximation of $f'(1)$ in the case $f(x)=e^x$ and $h=\frac{1}{2^i},\ i=1,2,\ldots,40$. Can you see the effect of round-off errors appearing for small values of $h$? Warning: by default mp.arange makes a vector of integers and python tries to compute powers of integers as integers, which gives incorrect results for large powers. Therefore, it may be a good idea to add an option dtype=float to np.arange command.
Often we can not compute the value of a function for all $x$ values, but have the values of the function available only at $n+1$ equally spaced points in an interval $[a,\,b]$. This means that we have a vector of values $(f_0,f_1,\ldots,f_n)$, where $f_i=f(x_i)$ for $x_i=a+i\cdot \frac{b-a}{n}$. In such case we still can use finite difference formulas for computing derivatives at the points $x_i,\ i=1,\ldots,n-1$, but we can not use arbitrary $h$ values. Instead, it makes sense to use the formulas with $h=\frac{b-a}{n}$. So, if we want to compute approximately $f'(c)$, where $c$ is a grid point inside $[a,\,b]$, we first have to find $i$ such that $c=x_i$ and then compute the approximate value of the derivative by $\frac{f_{i+1}-f_{i-1}}{2h}$.
Define a function myDer2 which for given vector f and for given numbers a,b,x computes approximately $f'(x)$. Assume that f corresponds to the uniform grid in the interval $[a,\,b]$ and that x corresponds to some grid point $x_i$. Check the correctness of your function in the case f=[1,0,1,4,9], a=-1, b=3 and x=1 (the answer should be 2).
def myDer2(f,a,b,x):
#compute the number of intervals n from the length of vector f
n=??
#compute h
h=???
#compute index i such that x[i]=x. Make sure that the result is an integer (use np.int64(), if needed)
i=???
return
myDer2(f=[1,0,1,4,9],a=-1,b=3,x=1)
Finite difference approximations of derivatives can be used for deriving numerical methods for solving differential equations.
Let us consider the following problem: find $y$ such that $$y''(x)+x\,y'(x)=f(x),\ x\in [a,b]$$ and $$y(a)=0,\ y(b)=1.$$ where $f$ is a given functions. The procedure for deriving a finite difference approximation for the problem above consists of the following steps.
If $x_i,\ i=0\ldots,n$ are equally spaced (with step size $h=\frac{b-a}{n}$), then we can use the finite difference approximation discussed above: $$y''(x_i)\approx \frac{y(x_{i-1})-2y(x_i)+y(x_{i+1})}{h^2}\qquad \text{(error}\leq ch^2),$$ $$y'(x_i)\approx \frac{y(x_{i+1})-y(x_{i-1})}{2h}\qquad \text{(error}\leq ch^2).$$ Hence, writing out the differential equation at points $x_i,\ i=1,\ldots,n-1$ and replacing derivatives with finite difference approximations we get a system of equations $$\frac{y_{i-1}-2y_i+y_{i+1}}{h^2}+x_i \frac{y_{i+1}-y_{i-1}}{2h}=f(x_i),\ i=1,\ldots,n-1,$$ where $y_i$ denotes the approximate values of $y(x_i)$. We can view the boundary conditions as two additional equations $$y_0=1,\ y_n=-1.$$ Since $i$-th equation contains only 3 unknowns $y_{i-1},y_i,y_{i+1}$ the resulting matrix of the system of equations has non-zero entries only on three diagonals, so we obtain a three diagonal system of equations.
It is natural to number the equations according to the grid points for which they were derived. So, the equation corresponding to the boundary condition at $x_0=a$ has number 0 (numbering starts from 0 in python!), the equation corresponding to $x_1$ has number 1 and so on. Let the coefficient of $y_j$ in the $i$-th equation be denoted by $m_{ij}$. These coefficient form the coefficient matrix of the system of equations we want to solve. Since the equation for $i=0$ is $y_0=1$, we have $m_{00}=1$ and $m_{0j}=0,\ j=1,2,\ldots,n$. For $i=1$ the equation is $$\frac{y_{0}-2y_1+y_{2}}{h^2}+x_1 \frac{y_{2}-y_{0}}{2h}=f(x_1)$$ and we get $m_{10}=\frac{1}{h^2}-\frac{x_1}{2h}$, $m_{11}=-\frac{2}{h^2}$, $m_{12}=\frac{1}{h^2}+\frac{x_1}{2h}$, $m_{1j}=0,\ j>2$. So it is quite easy to fill the coeficient matrix and it is quite easy to see that in the $i$-th equation only the coefficients $m_{i,i-1},m_{ii}$ and $m_{i,i+1}$ are different from 0 for $i=1,\ldots,n-1$.
Consider the following problem: find $y$ such that
$$y''(x)+xy'(x)=f(x),\ x\in [a,b]$$
$$y(a)=0,\ y(b)=1$$
One possibility to solve a linear system of equations is to use the Python command linalg.solve(M,z) from the subpackage linalg of SciPy (that has to be imported first with from scipy import linalg), which returns the solution $y$ of the system of equations $My=z$. Use this command to find the values of the approximate solution of the problem with the finite difference method in the case $n=10$, $a=0$, $b=1$, $f(x)=3x^3+6x$. Find the errors between the approximate solution and the exact solution $y(x)=x^3$.
Solution
#your code here!
If the system matrix has only some nonzero diagonals then it is actually a waste of computer memory to store the full matrix. For solving such system it is actually possible to use the command linalg.solve_banded((l,u),Dgs,z), where the rows of the matrix $Dgs$ contain the diagonals of $M$ starting from the highest one, $l$ is the number of diagonals below the main diagonal and $u$ is the number of diagonals above the main diagonal (in our case $l=u=1$).NB! In the matrix Dgs for the diagonals that are above the main diagonal the first elements are actually not used (for the diagonal directly above the main diagonal the first element is not used, for the next diagonal two steps above the main diagonal the first two elements are not used etc); for diagonals below the main diagonal last elements are not used. Write a function that for a given $n$ solves the problem of the previous exercise with the command linalg.solve_banded((l,u),Dgs,z) and returns the maximal error at the points $x_i,\ i=1,\ldots, n-1$. Determine how many times the error is reduced if we increase the number of points two times.
Solution
Practical homework 3 is related to derivation of a finite difference method for solving a boundary value problem of an ordinary differential equation (similar to lab exercise 4). The Moodle VPL exercise "Practical homework 3" is used for assigning problems and for checking solutions. Each student is randomly assigned a variant of the equation.