Numerical Calculus
Contents
5. Numerical Calculus#
Throughout this section and the next ones, we shall cover the topic of numerical calculus. Calculus has been identified since ancient times as a powerful toolkit for analysing and handling geometrical problems. Since differential calculus was developed by Newton and Leibniz (in its actual notation), many different applications have been found, at the point that most of the current science is founded on it (e.g. differential and integral equations). Due to the ever increasing complexity of analytical expressions used in physics and astronomy, their usage becomes more and more impractical, and numerical approaches are more than necessary when one wants to go deeper. This issue has been identified since long ago and many numerical techniques have been developed. We shall cover only the most basic schemes, but also providing a basis for more formal approaches.
5.1. Bibliography#
https://github.com/restrepo/Calculus
Thomas J. Sargent John Stachurski, Python Programming for Quantitative Economics GitHub
Numerical Integration
%pylab inline
import numpy as np
import scipy.interpolate as interp
%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib
5.2. Numerical Differentiation#
According to the formal definition of differentiation, given a function \(f(x)\) such that \(f(x)\in C^1[a,b]\), the first order derivative is given by
However, when \(f(x)\) exhibits a complex form or is a numerical function (only a discrete set of points are known), this expression becomes unfeasible. In spite of this, this formula gives us a very first rough way to calculate numerical derivatives by taking a finite interval \(h\), i.e.
where the function must be known at least in \(x_0\) and \(x_1 = x_0+h\), and \(h\) should be small enough.
5.2.1. Example 1#
Evaluate the first derivative of the next function using the previous numerical scheme at the point \(x_0=2.0\) and using \(h=0.5,\ 0.1,\ 0.05,...\)
\(f(x) = \sqrt{1+\cos^2(x)}\)
Compare with the real function and plot the tangent line using the found values of the slope.
from scipy import misc
It is used as:
misc.derivative(func, x0, dx=1.0, n=1, args=(), order=3)
Parameters
func
: function →
Input function.
x0
: float →
The point at which n
-th derivative is found.
dx
: float, optional →
Spacing.
n
: int, optional →
Order of the derivative. Default is 1.
args
: tuple, optional →
Arguments
order
: int, optional →
Number of points to use, must be odd.
def f(x):
return np.sqrt( 1+np.cos(x)**2 )
x0 = 2
Activity: We now check for the impact of the change of the spacing dx
. Try from dx=0.5
and then small values
eval('np.cos(np.pi)/np.sqrt(2)')
-0.7071067811865475
for i in range(5):
dx=eval(input('dx='))
print( misc.derivative(f,x0,dx=dx) )
dx= 1
0.13526274947015327
dx= 0.1
0.3462499420237386
dx= 0.01
0.3493266965195363
dx= 1e-3
0.34935758881293744
dx= 1e-6
0.3493579009417047
Compare with:
misc.derivative(f,x0,dx=1E-8)
0.34935789816614715
fp=lambda x:-np.cos(x)*np.sin(x)/f(x)
fp(2)
0.3493579008690488
m=misc.derivative(f,x0,dx=1E-6)
b=f(x0)-m*x0
xmin = 1.6
xmax = 2.4
x = np.linspace( xmin, xmax, 100 )
plt.plot( x, f(x), color="black", label="function", linewidth=1 )
plt.plot( x, m*x+b,"r--", label=f"derivative at $x={x0}$")
plt.plot(x0,f(x0),'y*',markersize=10)
plt.grid()
plt.legend(fontsize=15)
plt.xlabel('$x$',size=15)
plt.ylabel(r'$f(x) = \sqrt{1+\cos^2(x)}$',size=15)
Text(0, 0.5, '$f(x) = \\sqrt{1+\\cos^2(x)}$')
5.2.2. Implementation of the derivate of the function inside a full range#
We now generalize the derivative
function to allow the evaluation of the derivate in a full range of values. It will be designed such that the evaluation in just one point can be still possible. In this way, the new function can be used as a full replacement of derivative
function.
To implement this function we need three keys ingredients of python
:
try
andexcept
python progamming sctructureFunction definition with generic mandatory and optional arguments
Conversion of a float function into a vectorized fuction
5.2.2.1. try
and except
python progamming sctructure#
First we introduce the try
and except
python progamming sctructure, wich is used to bypass one python error. For example a zero dimension array has not a shape attribute, so that the following error, of type IndexError
, is generated:
(1,)[1]
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
Input In [20], in <cell line: 1>()
----> 1 (1,)[1]
IndexError: tuple index out of range
np.array(3)
array(3)
nn=np.array(3).shape[0]
nn
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
Input In [33], in <cell line: 1>()
----> 1 nn=np.array(3).shape[0]
2 nn
IndexError: tuple index out of range
To bypass that error we use the following code:
x=3
len(3)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
/tmp/ipykernel_57720/4138546388.py in <module>
1 x=3
----> 2 len(3)
TypeError: object of type 'int' has no len()
np.array([]).shape
(0,)
x='3' #[2,3] #3 #[2,3], array([])
try:
nn=np.array(x).shape[0]
except IndexError:
nn=-1
so that nn
takes the values assigned in the except
part:
nn
2
Activity: A float has the method is_integer()
to check if the decimal part is zero. Creates a function is_integer(x)
which works even in the case that the number is already an integer
x=2.
x.is_integer()
True
x=2
x.is_integer()
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Input In [50], in <cell line: 2>()
1 x=2
----> 2 x.is_integer()
AttributeError: 'int' object has no attribute 'is_integer'
x=2.3
x.is_integer()
False
isinstance(2,int)
True
Objetos:
str,int,float,list,dict
isinstance(2,list)
False
eval('hola')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [58], in <cell line: 1>()
----> 1 eval('hola')
File <string>:1, in <module>
NameError: name 'hola' is not defined
eval('hola mundo')
Traceback (most recent call last):
File ~/anaconda3/lib/python3.9/site-packages/IPython/core/interactiveshell.py:3369 in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
Input In [59] in <cell line: 1>
eval('hola mundo')
File <string>:1
hola mundo
^
SyntaxError: unexpected EOF while parsing
def is_integer(x):
'''
Returns true if x is integer even if x is integer
'''
try:
return x.is_integer()
except AttributeError:
if isinstance(x,int):
return True
elif isinstance(x,str):
try:
#recursive
return is_integer(eval(x))
except (NameError,SyntaxError):
return False
else:
return False
## Test the functions for typical cases
assert is_integer(3.)
assert is_integer(3.5)==False
assert is_integer(3)
assert is_integer('3')
assert is_integer('3+9+11')
assert is_integer('3.5')==False
assert is_integer('a')==False
assert is_integer('a b')==False
assert is_integer('1/2+1/2')
assert is_integer([1,3,5])==False
assert is_integer({'A':1})==False
5.2.3. Testing#
assert True
assert False
---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
Input In [61], in <cell line: 1>()
----> 1 assert False
AssertionError:
assert is_integer(3)
True
assert
is trivial for True
assert
generates an error for False
5.2.3.1. Function parameters#
In python the function can be defined with mandatory and optional arguments
def f(a,b,c=2,d=3):
return a+b-c+d
a
,b
: is like a tuple of mandatory arguments →(a,b)
c=2
,d=3
: is like a dictionary of optional arguments →{'c':2,'d':3}
5.2.3.2. Function definition with generic mandatory and optional arguments#
We need to introduce another important concept about functions. There is a powerfull way to define a function with generic arguments:
Instead of the mandatory arguments we can use the generic list pointer:
*args
Instead of the optional arguments we can use the generic dictionary pointer:
**kwargs
def f(*args,**kwargs):
if args:
print('*args is the tuple of mandatory arguments: {}'.format(args))
print('This allows to have a dynamic return according to input:')
return type(args[0])
if kwargs:
print('*kwargs is the dictionary of optional arguments: {}'.format(kwargs))
print('This allows to have a dynamic return according to input:')
return type(list(kwargs.keys())[0]),type(list(kwargs.values())[0])
Activity; Check the function with any type of mandatory or optional argument
Check the function without arguments
f()
Check the function with a mandatory argument
f('hola mundo')
*args is the tuple of mandatory arguments: ('hola mundo',)
This allows to have a dynamic return according to input:
str
Check the function with several mandatory arguments
f(2,[1,2],1+3j)
*args is the tuple of mandatory arguments: (2, [1, 2], (1+3j))
This allows to have a dynamic return according to input:
int
Check the function with several optional arguments
Check the function with one of the optional arguments equal to some list
f(d=3,h=[1,2],k='kk')
*kwargs is the dictionary of optional arguments: {'d': 3, 'h': [1, 2], 'k': 'kk'}
This allows to have a dynamic return according to input:
(str, int)
Check the function with one string as mandatory argument, and with one of the optional arguments equal to some list
5.2.3.3. Conversion of a float function into a vectorized function#
An important numpy
function is:
Which convert a function which takes a float as an argument into a new function which que takes numpy arrays
as an argument. The problem is that the converted function does not longer return a float when the argument is a float:
It is used as:
np.vectorize(pyfunc,...)
Parameters:
pyfunc
: A python function or method.
…
import math as m
m.sin(0.5)
0.479425538604203
m.sin([0.5,0.7])
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Input In [80], in <cell line: 1>()
----> 1 m.sin([0.5,0.7])
TypeError: must be real number, not list
import math as m
sin=np.vectorize(m.sin)
sin([0.5,0.7])
array([0.47942554, 0.64421769])
sin(0.5)
array(0.47942554)
5.2.3.4. Vectorization of Scipy method derivative
#
The problem is that derivative
only works for one a point
misc.derivative(np.sin,1,dx=1E-6)
0.5403023058958567
misc.derivative(np.sin,[1,1.2],dx=1E-6)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Input In [84], in <cell line: 1>()
----> 1 misc.derivative(np.sin,[1,1.2],dx=1E-6)
File ~/anaconda3/lib/python3.9/site-packages/scipy/misc/_common.py:144, in derivative(func, x0, dx, n, args, order)
142 ho = order >> 1
143 for k in range(order):
--> 144 val += weights[k]*func(x0+(k-ho)*dx,*args)
145 return val / prod((dx,)*n,axis=0)
TypeError: can only concatenate list (not "float") to list
This is easily fixed with np.vectorize
derivative=np.vectorize(misc.derivative)
derivative(np.sin,[1,1.2],dx=1E-6)
array([0.54030231, 0.36235775])
derivative(np.sin,1,dx=1E-6)
array(0.54030231)
misc.derivative(np.sin,1,dx=1E-6)
0.5403023058958567
5.2.4. Code for the derivate of the function inside a full range#
To recover the evaluation of a float into a float we can force the IndexError
in order to use of the pure scipy
derivative
function when a float is given as input. The full implemention combine the three previos ingredients into a very compact pythonic-way function definition as shwon below
%pylab inline
import numpy as np
from scipy import misc
def derivative(func,x0,**kwargs):
'''
Vectorized replacement of scipy.misc derivative:
from scipy.misc import derivative
For usage check the derivative help, e.g, in jupyter:
from scipy.misc import derivative
derivative?
'''
try:
#x0: can be an array or a list
nn=np.asarray(x0).shape[0] ## force error if float is used
fp=np.vectorize(misc.derivative)
except IndexError:
fp=misc.derivative
return fp(func,x0,**kwargs)
assert isinstance(derivative(np.sin,1,dx=1E-6),float)
import numpy as np
A=np.array([2])
type(A[0])
numpy.int64
numbers.Real
numbers.Real
import numbers
x=2#A[0]
isinstance(x, numbers.Real)#numbers.Integral)
True
import numpy as np
from scipy import misc
def derivative(func,x0,**kwargs):
'''
Vectorized replacement of scipy.misc derivative:
from scipy.misc import derivative
For usage check the derivative help, e.g, in jupyter:
from scipy.misc import derivative
derivative?
'''
try:
#x0: can be an array or a list
nn=np.asarray(x0).shape[0] ## force error if float is used
fp=np.vectorize(misc.derivative)
except IndexError:
fp=misc.derivative
return fp(func,x0,**kwargs)
assert isinstance(derivative(np.sin,1,dx=1E-6),float)
which behaves exactly as the scipy
derivative
function when a float argument is used:
derivative(np.sin,1,dx=1E-6)
0.5403023058958567
derivative(np.sin,[1,1.5],dx=1E-6)
array([0.54030231, 0.0707372 ])
and as a vectorized function when an array argument is used:
misc.derivative(np.sin,1,dx=1E-6),misc.derivative(np.sin,1.5,dx=1E-6)
(0.5403023058958567, 0.0707372017072494)
derivative(np.sin,[1.,1.5],dx=1E-6)
array([0.54030231, 0.0707372 ])
derivative(np.sin,[1.,1.5],dx=1E-6,n=4,order=5)
array([-7.77156117e+08, -1.11022302e+09])
Let see now the implementation of the derivate
function as full replacement of the derivative
function:
5.2.4.1. Test#
Let us check the implementation with the previous function
func=lambda x: np.sqrt( 1+np.cos(x)**2 )
but now evaluated for a list of values
x=[1.8,2,2.2]
funcp=derivative(func,x,dx=1E-3)
funcp
array([0.21576117, 0.34935759, 0.41006125])
For the prevous X
array:
x0 = 2
xmin = 1.6
xmax = 2.4
x = np.linspace( xmin, xmax, 100 )
We have:
derivative(func,x,dx=1E-3).shape
(100,)
xmin = 1.6
xmax = 2.4
X = np.linspace( xmin, xmax, 100 )
plt.plot( X, func(X) , color="black", label="$f(x)$", linewidth=3)
plt.plot( X, derivative(func,X,dx=1E-3), color="red", label="$f'(x)$", linewidth=3)
plt.plot( X, derivative(func,X,dx=1E-3,n=2), color="blue", label="$f''(x)$", linewidth=3)
plt.legend(loc='best',fontsize=15)
plt.xlabel('$x$',size=15)
plt.grid()
5.2.4.2. Example#
Finally we check the function \(f(x)=\cos x\) in the range \([0,2\pi]\)
xmin = 0
xmax = 2*np.pi
X = np.linspace( xmin, xmax, 100 )
plt.plot( X, np.cos(X) , "g-", label="$f(x)=\cos x$", linewidth=3)
plt.plot( X, derivative(np.cos,X,dx=1E-3), color="red", label="$f'(x)$", linewidth=3)
plt.plot( X, -np.sin(X) , 'b:', label="$-\sin x$", linewidth=3, zorder=10 )
plt.plot( X, derivative(np.cos,X,dx=1E-3,n=4,order=5), 'k:', label="$f^{(iv)}(x)$", linewidth=3)
plt.legend(loc='best',fontsize=15)
plt.xlabel('$x$',size=15)
Text(0.5, 0, '$x$')
misc.derivative(np.cos,np.pi/2,n=4,dx=1E-3,order=5)
-1.3010426069826051e-06
Activity: Implement the full derivative function by using isinstance()
instead of try
and except
from scipy import misc
def derivate(func,x0,**kwargs):
'''
Vectorized replacement of scipy.misc derivative:
from scipy.misc import derivative
For usage check the derivative help, e.g, in jupyter:
from scipy.misc import derivative
derivative?
'''
if isinstance(x0,list): ## or ... or ...
print(x0)
fp=np.vectorize(misc.derivative)
else:
fp=misc.derivative
return fp(func,x0,**kwargs)
isinstance( np.array([2,3]),numpy.ndarray )
True
derivate(np.cos,np.array([2,3]),dx=1E-6)
array([-0.90929743, -0.14112001])
5.2.5. Example: Heat transfer in a 1D bar#
Fourier’s Law of thermal conduction describes the diffusion of heat. Situations in which there are gradients of heat, a flux that tends to homogenise the temperature arises as a consequence of collisions of particles within a body. The Fourier’s Law is giving by
where T is the temperature, \(\nabla T\) its gradient and k is the material’s conductivity. In the next example it is shown the magnitud of the heat flux in a 1D bar(wire).
def Temp(x):
return x**3 + 3*x-1
Xn = np.linspace(0,10,100)
#Temperature profile
def Temp(x):
return x**3 + 3*x-1
## Points where function is known
Xn = np.linspace(0,10,100)
plt.figure( figsize=(8,7) )
plt.plot(Xn,derivate(Temp,Xn)*10)
plt.grid()
plt.xlabel( "$x$",fontsize =15 )
plt.ylabel( r"$\frac{dT}{dx}$",fontsize =20 )
plt.title( " Magnitud heat flux transfer in 1D bar" )
Text(0.5,1,' Magnitud heat flux transfer in 1D bar')
5.2.6. Activity#
Construct a density map of the magnitud of the heat flux of a 2D bar. Consider the temperature profile as $\( T(x,y) = x^3 + 3x-1+y^2 \)$
5.2.7. Activity#
The Poisson’s equation relates the matter content of a body with the gravitational potential through the next equation
where \(\phi\) is the potential, \(\rho\) the density and \(G\) the gravitational constant.
Taking these data and using the three-point Midpoint formula, find the density field from the potential (seventh column in the file) and plot it against the radial coordinate. (Tip: Use \(G=1\))
5.2.8. Activity#
The radar stations A and B, separated by the distance a = 500 m, track the plane C by recording the angles \(\alpha\) and \(\beta\) at 1-second intervals. The successive readings are
calculate the speed v using the 3 point approximantion at t = 10 ,12 and 14 s. Calculate the x component of the acceleration of the plane at = 12 s. The coordinates of the plane can be shown to be
5.3. Numerical Integration#
5.4. Appendix#
5.4.1. One implementation of derivative
algorithm#
#Function to evaluate
def function(x):
return np.sqrt( 1+np.cos(x)**2 )
#X value
x0 = 2
xmin = 1.8
xmax = 2.2
#h step
hs = [0.5,0.1,0.05]
#Calculating derivatives
dfs = []
for h in hs:
dfs.append( (function(x0+h)-function(x0))/h )
#Plotting
plt.figure( figsize=(10,8) )
#X array
X = np.linspace( xmin, xmax, 100 )
Y = function(X)
plt.plot( X, Y, color="black", label="function", linewidth=3, zorder=10 )
#Slopes
Xslp = [1,x0,3]
Yslp = [0,0,0]
for df, h in zip(dfs, hs):
#First point
Yslp[0] = function(x0)+df*(Xslp[0]-Xslp[1])
#Second point
Yslp[1] = function(x0)
#Third point
Yslp[2] = function(x0)+df*(Xslp[2]-Xslp[1])
#Plotting this slope
plt.plot( Xslp, Yslp, linewidth = 2, label="slope$=%1.2f$ for $h=%1.2f$"%(df,h) )
#Format
plt.grid()
plt.xlabel("x")
plt.ylabel("y")
plt.xlim( xmin, xmax )
plt.ylim( 1, 1.15 )
plt.legend( loc = "upper left" )
<matplotlib.legend.Legend at 0x7f2790fe7990>
5.4.2. n+1-point formula#
A generalization of the previous formula is given by the (n+1)-point formula, where first-order derivatives are calculated using more than one point, what makes it a much better approximation for many problems. it is controled by the order
option of the derivative
function of scipy.misc
Theorem
For a function \(f(x)\) such that \(f(x)\in C^{n+1}[a,b]\), the next expression is always satisfied
where \(\{x_i\}_i\) is a set of point where the function is mapped, \(\xi(x)\) is some function of \(x\) such that \(\xi\in[a,b]\), and \(P(x)\) is the associated Lagrange interpolant polynomial.
As \(n\) becomes higher, the approximation should be better as the error term becomes neglectable.
Taking the previous expression, and differenciating, we obtain
where \(L_{n,k}\) is the \(k\)-th Lagrange basis functions for \(n\) points, \(L'_{n,k}\) is its first derivative.
Note that the last expressions is evaluated in \(x_j\) rather than a general \(x\) value, the cause of this is because this expression is not longer valid for another value not within the set \(\{x_i\}_i\), however this is not an inconvenient when handling real applications.
This formula constitutes the (n+1)-point approximation and it comprises a generalization of almost all the existing schemes to differentiate numerically. Next, we shall derive some very used formulas.
For example, the form that takes this derivative polynomial for 3 points \((x_i,y_i)\) is the following
5.4.3. Endpoint formulas#
Endpoint formulas are based on evaluating the derivative at the first of a set of points, i.e., if we want to evaluate \(f'(x)\) at \(x_i\), we then need \((x_i\), \(x_{i+1}=x_i+h\), \(x_{i+2}=x_i+2h\), \(\cdots)\). For the sake of simplicity, it is usually assumed that the set \(\{x_i\}_i\) is equally spaced such that \(x_k = x_0+k\cdot h\).
Three-point Endpoint Formula
with \(\xi\in[x_i,x_i+2h]\)
Five-point Endpoint Formula
with \(\xi\in[x_i,x_i+4h]\)
Endpoint formulas are especially useful near to the end of a set of points, where no further points exist.
5.4.4. Midpoint formulas#
On the other hand, Midpoint formulas are based on evaluating the derivative at the middle of a set of points, i.e., if we want to evaluate \(f'(x)\) at \(x_i\), we then need \((\cdots\), \(x_{i-2} = x_i - 2h\), \(x_{i-1} = x_i - h\), \(x_i\), \(x_{i+1}=x_i+h\), \(x_{i+2}=x_i+2h\), \(\cdots)\).
Three-point Midpoint Formula
with \(\xi\in[x_i-h,x_i+h]\)
Five-point Midpoint Formula
with \(\xi\in[x_i-2h,x_i+2h]\)
As Midpoint formulas required one iteration less than Endpoint ones, they are more often used for numerical applications. Furthermore, the round-off error is smaller as well. However, near to the end of a set of points, they are no longer useful as no further points exists, and Endpoint formulas are preferable.
Example with custom implementation
#Temperature profile
def Temp(x):
return x**3 + 3*x-1
#Derivative three end point
def TEP( Yn,i, h=0.01,right=0 ):
suma = -3*Yn[i]+4*Yn[i+(-1)**right*1]-Yn[i+(-1)**right*2]
return suma/(2*h*(-1)**right)
#Derivative mid point
def TMP( Ynh,Ynmh, h = 0.01 ):
return (Ynh-Ynmh)/(2*h)
## Points where function is known
Xn = np.linspace(0,10,100)
Tn = Temp(Xn)
#Magnitude of heat flux array
Q = np.zeros(len(Xn))
#Left end derivative
Q[0] = TEP(Tn,0)
#Mid point derivatives
index = len(Xn)-1
for i in xrange( 1,index ):
Q[i] = TMP( Tn[i+1],Tn[i-1] )
#Right end derivative
Q[-1] = TEP( Tn,index,right=1 )
#Plotting
plt.figure( figsize=(8,7) )
plt.plot(Xn,Q)
plt.grid()
plt.xlabel( "x",fontsize =15 )
plt.ylabel( "$\\frac{dT}{dx}$",fontsize =20 )
plt.title( " Magnitud heat flux transfer in 1D bar" )
<matplotlib.text.Text at 0x7f87263de3d0>