Python for Engineering
Tim Pengajar IF2132 Sem. 1 2018/2019
© 2018
Why Python for Scientific Computing?
The design focus on the Python language is on productivity and code
readability, for example through:
• Interactive python console
• Very clear, readable syntax through whitespace indentation
• Strong introspection capabilities
• Full modularity, supporting hierarchical packages
• Exception-based error handling
• Dynamic data types & automatic memory management
10/31/2018 IF2132/Pemprograman Komputer 2
Why Python for Scientific Computing?
•Optimisation strategies
•Get it right first, then make it fast
•Prototyping in Python
10/31/2018 IF2132/Pemprograman Komputer 3
Numerical Python
The NumPy package (read as NUMerical PYthon) provides access to
• a new data structure called arrays which allow
• efficient vector and matrix operations. It also provides
• a number of linear algebra operations (such as solving of systems of
linear equations, computation of Eigenvectors and Eigenvalues).
10/31/2018 IF2132/Pemprograman Komputer 4
Vector
>>> import numpy as N >>> c = N. zeros (4)
>>> a = N. array ([0 , >>> c [0] = 3.4
0.5 , 1, 1.5]) >>> c [2] = 4
>>> print (a) >>> print (c)
>>> print (c [0])
>>> b = N. arange (0, 2, >>> print (c [0: -1])
0.5)
>>> print (b)
10/31/2018 IF2132/Pemprograman Komputer 5
Matrix Multiplication
>>> import numpy as N
>>> import numpy . random
>>> A = numpy . random . rand (5, 5) # generates a random 5 by 5 matrix
>>> x = numpy . random . rand (5) # generates a 5- element vector
>>> b=N. dot (A, x) # multiply matrix A with vector x
10/31/2018 IF2132/Pemprograman Komputer 6
Solving systems of linear equations
>>> import numpy . linalg as LA
>>> x = LA. solve (A, b)
10/31/2018 IF2132/Pemprograman Komputer 7
Computing Eigenvectors and Eigenvalues
>>> import numpy
>>> import numpy . linalg as LA
>>> A = numpy .eye (3) #'eye '->I ->1 (ones on the diagonal)
>>> print (A)
>>> evalues , evectors = [Link] (A)
>>> print ( evalues )
>>> print ( evectors )
10/31/2018 IF2132/Pemprograman Komputer 8
Curve fitting of polynomials
# demo curve fitting : xdata and ydata are input data # evaluate p(x) for all x in list xs
xdata = numpy . array ([0.0 , 1.0 , 2.0 , 3.0 , 4.0 , 5.0]) import pylab
ydata = numpy . array ([0.0 , 0.8 , 0.9 , 0.1 , -0.8, -1.0]) pylab . plot (xdata , ydata , 'o', label ='data')
z = numpy . polyfit (xdata , ydata , 3) pylab . plot (xs , ys , label ='fitted curve')
p = numpy . poly1d (z) pylab . ylabel ('y')
pylab . xlabel ('x')
# create plot pylab . savefig ('[Link] ')
pylab . show ()
xs = [0.1 * i for i in range (50)]
ys = [p(x) for x in xs]
10/31/2018 IF2132/Pemprograman Komputer 9
Fine tunning your plot
import pylab pylab . plot (x, y2 , label ='cos(x)')
pylab . legend ()
import numpy as N
pylab . grid ()
x = N. arange ( -3.14 , 3.14 , 0.01)
pylab . xlabel ('x')
y1 = N. sin(x) pylab . title ('This is the title of the graph')
y2 = N. cos(x) pylab . show () # necessary to display graph
pylab . figure ( figsize =(5 , 5))
pylab . plot (x, y1 , label ='sin(x)')
10/31/2018 IF2132/Pemprograman Komputer 10
Numerical integration
from numpy import *
from [Link] import quad
# function we want to integrate
def f(x):
return x*x
# call quad to integrate f from -2 to 2
res , err = quad (f, 1, 2)
print ("The numerical result is {:f} (+ -{:g})".format (res , err ))
10/31/2018 IF2132/Pemprograman Komputer 11
Solving Ordinary Differential Equation
from scipy . integrate import odeint y = odeint (f, y0 , t)
import numpy as N # actual computation of y(t)
def f(y, t):
import pylab # plotting of results
return -2 * y * t
pylab . plot (t, y)
y0 = 1 # initial value
pylab . xlabel ('t'); pylab . ylabel ('y(t)')
a = 0 # integration limits for t
pylab . show ()
b = 2
t = N. arange (a, b, 0.01)
# values of t for which we require
# the solution y(t)
10/31/2018 IF2132/Pemprograman Komputer 12
Root Finding using Bisection Method
from scipy . optimize import bisect
def f(x):
""" returns f(x)=x^3 -2x^2. Has roots at
x=0 ( double root ) and x=2 """
return x ** 3 - 2 * x ** 2
# main program starts here
x = bisect (f, 1.5 , 3, xtol =1e-6)
print ("The root x is approximately x =%14.12g ,\n"
"the error is less than 1e -6." % (x))
print ("The exact error is %g." % (2 - x))
10/31/2018 IF2132/Pemprograman Komputer 13
Interpolation
# main program
n = 10
x, y = create_data (n)
import numpy as np #use finer and regular mesh for plot
xfine = np. linspace (0.1 , 4.9 , n * 100)
import scipy . interpolate
# interpolate with piecewise constant function (p=0)
import pylab
y0 = scipy . interpolate . interp1d (x, y, kind ='nearest')
def create_data (n): # interpolate with piecewise linear func (p=1)
xmax = 5. y1 = scipy . interpolate . interp1d (x, y, kind ='linear')
x = np. linspace (0, xmax , n) # interpolate with piecewise constant func (p=2)
y = - x**2 y2 = scipy . interpolate . interp1d (x, y, kind ='quadratic')
y += 1.5 * np. random . normal ( size =len(x))
return x, y
10/31/2018 IF2132/Pemprograman Komputer 14
Interpolation (2)
pylab . plot (x, y, 'o', label ='data point')
pylab . plot (xfine , y0( xfine ), label ='nearest')
pylab . plot (xfine , y1( xfine ), label ='linear')
pylab . plot (xfine , y2( xfine ), label ='cubic')
pylab . legend ()
pylab . xlabel ('x')
pylab . savefig (' [Link]')
pylab . show ()
10/31/2018 IF2132/Pemprograman Komputer 15
Curve Fitting
import numpy as np # call curve fit function
popt , pcov = curve_fit (f, x, yi)
from scipy . optimize import curve_fit
a, b, c = popt
def f(x, a, b, c): print (" Optimal parameters are a=%g, b=%g, and c=%g" % (a, b, c))
""" Fit function y=f(x,p) with parameters # plotting
p=(a,b,c). """ import pylab
return a * np. exp (- b * x) + c yfitted = f(x, * popt ) # equivalent to f(x, popt [0] , popt [1] , popt
[2])
# create fake data pylab . plot (x, yi , 'o', label ='data $y_i$ ')
x = np. linspace (0, 4, 50) pylab . plot (x, yfitted , '-', label ='fit $f(x_i)$')
y = f(x, a=2.5 , b=1.3 , c =0.5) pylab . xlabel ('x')
pylab . legend ()
#add noise
pylab . savefig ('[Link]')
yi = y + 0.2 * np. random . normal ( size =len (x))
10/31/2018 IF2132/Pemprograman Komputer 16
