Python for Numerical Computing and Scientific Applications
P
r
o
f
.
F
e
d
e
r
i
c
o
C
l
u
n
i
UNIVERSITÀ DEGLI STUDI DI PERUGIA
International Doctoral Program in Civil And Environmental Engineering
A
.
Y
.
2
0
2
0
/
2
1
L
e
s
s
o
n
#
3
–
S
e
p
t
e
m
b
e
r
,
2
1
t
h
2
0
2
1
B
A
S
E
I
N
T
E
R
A
C
T
I
V
I
T
Y
C
o
m
m
a
n
d
l
i
n
e
i
n
t
e
r
f
a
c
e
i
n
t
e
r
a
c
t
i
v
i
t
y
It is possible to build a script which interact with the user asking for specific data
def radici(a,b,c):
delta = b**2-4*a*c
x1 = (-b-delta**0.5)/(2*a)
x2 = (-b+delta**0.5)/(2*a)
return x1, x2
a = float(input('Inserisci a: '))b = float(input('Inserisci b: '))c = float(input('Inserisci c: '))x1, x2 = radici(a,b,c)
print('x1 = '+str(x1))print('x2 = '+str(x2))N.B “input” read a
string
from
the command line
B
A
S
E
I
N
T
E
R
A
C
T
I
V
I
T
Y
C
o
m
m
a
n
d
l
i
n
e
i
n
t
e
r
f
a
c
e
i
n
t
e
r
a
c
t
i
v
i
t
y
Moreover, it is also possible to execute a script from the OS shell, and eventually to
pass arguments directly from the shell
import sys
def radici(a,b,c):
delta = b**2-4*a*c
x1 = (-b-delta**0.5)/(2*a)
x2 = (-b+delta**0.5)/(2*a)
return x1, x2
if __name__=='__main__':
if len(sys.argv)>1:
a = float(sys.argv[1])
b = float(sys.argv[2])
c = float(sys.argv[3])
else:
a = float(input('Inserisci a: ')) b = float(input('Inserisci b: ')) c = float(input('Inserisci c: '))x1, x2 = radici(a,b,c)
print('x1 = '+str(x1)) print('x2 = '+str(x2))radici.py
> python radici.py
> python radici.py 1 1 -2
T
H
E
S
C
I
P
Y
“
E
C
O
S
Y
S
T
E
M
”
N
U
M
P
Y
A
R
R
A
Y
<
t
y
p
e
‘
n
d
a
r
r
a
y
’
>
It is a multi-dimensional (may have more than one index) with homogeneous
elements and fixed dimension.
There are constructors for some typical cases:
import numpy as np
a
= np.zeros(4)
a
= np.ones((2,3,4))
a
= np.arange(-5.,5.,1.)
a
= np.linspace(-5.,5.,100)
Array can also be generated through list and tuples:
l
= [1., 2., 3. ]
a
= np.array(
l
)
R
e
m
a
r
k
As usual with object,
b=a
create an alias and not a copy. To obtain a copy:
b
=
a
.copy()
N.B this is the “standard” way
of importing numpy
N
U
M
P
Y
A
R
R
A
Y
<
t
y
p
e
‘
n
d
a
r
r
a
y
’
>
It is possible to redistribute the elements among dimension.
import numpy.random as rd
a = rd.randint(-10,10,(2,3,4))
a.shape = (2,12)
or
a.re
shape (
2,12
)
Usual slicing operation to extract portions of the array:
a
[
indexes dim. 0
,
indexes dim. 1, :
]
N
o
t
e
t
h
a
t
s
l
i
c
i
n
g
p
r
o
d
u
c
e
v
i
e
w
s
o
f
t
h
e
a
r
r
a
y
,
i
.
e
.
a
[
:
,
:
]
d
o
e
s
n
o
t
p
r
o
d
u
c
e
a
n
e
w
a
r
r
a
y
!
T
o
t
a
k
e
a
l
l
e
l
e
m
e
n
t
s
o
f
o
n
e
d
i
m
e
n
s
i
o
n
,
u
s
e
«
:
»
To access the element at place i
a.item(i)
The 1-dimensional representation (as a view) is given by
a.ravel()
N.B numpy.random module is
more powerful than random
A
R
R
A
Y
<
t
y
p
e
‘
n
d
a
r
r
a
y
’
>
To iterate on the dimension
for i in range(a.shape[0]):
for j in range(a.shape[1]):
...
The usual operations
+
,
-
,
/
,
*
,
**
are possible, in such a case the operation is applied to all elements in the array
b = 2*a
Some useful methods
a.max() a.argmax()
a.min() a.argmin()
a.mean() a.std()
N
U
M
P
Y
M
A
T
R
I
X
<
t
y
p
e
‘
m
a
t
r
i
x
’
>
To define matrices:
m
= matrix(
list/array
)
or:
m
= mat(
list/array
)
Usual matrix operations are defined (in particular,
*
is matrix multiplication).
To perform matrix multiplication using arrays, use
@
.
In old Python versions, to matrix multiply two array
a
and
b
:
c=dot(a,b)
.
N
U
M
P
Y
T
h
e
“
R
o
s
e
t
t
a
s
t
o
n
e
”
f
o
r
M
a
t
l
a
b
u
s
e
r
s
https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
N
U
M
P
Y
An array can be used as the argument of a function. Anyway, if there is a
if
condition, the condition can not be evaluated:
def funz(x):
if x < 0:
return -1
else:
return 2*x
In such cases, the command
where
can be used:
def funz(x):
return np.where( x < 0, -1, 2*x)
It is worth noting that this last function can be called with a scalar as an argument.
Otherwise
funzv = np.vectorize(funz)
N
U
M
P
Y
-
V
E
C
T
O
R
I
A
L
I
Z
A
T
I
O
N
N
U
M
P
Y
-
U
S
E
F
U
L
S
U
B
M
O
D
U
L
E
S
/
F
U
N
C
T
I
O
N
S
There are several submodules in Numpy, for example random, fft, linalg
To import a submodule
import numpy.random as rd
import numpy.linalg as la
A = rd.randint(-10,10,(3,3))
b = rd.randint(-10,10,3)
if not np.isclose(la.det(A),0,atol=1e-12):
c = la.solve(A,b)
c1 = np.dot(la.inv(A),b)
w, v = la.eig(
A
) # solution of A@v = wv
N
U
M
P
Y
-
U
S
E
F
U
L
S
U
B
M
O
D
U
L
E
S
/
F
U
N
C
T
I
O
N
S
To write/read to a file (in text format)
import numpy as np
x = np.random.uniform(size=(10,3))
np.savetxt('dati.txt',x,'%7.4f',';')y = np.loadtxt('dati.txt’,delimiter=';')P
Y
L
A
B
/
M
A
T
P
L
O
T
L
I
B
To create graphs MatLab-like.
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
x = np.linspace(0.,2*np.pi, 101)
y = np.sin(x)
ys = x-x**3/6+x**5/120-x**7/5040+x**9/362880
plt.plot( x, y, 'bo-', label='true')
plt.plot( x, ys, 'r-', label='serie 5 terms')
plt.xlabel('$t$')plt.ylabel('$\sin(t^1)$')plt.title('Graph of $\sin(t)$')plt.legend()
plt.show()
To see several examples of use
https://matplotlib.org/gallery/index.html
To plot 3d graphs
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
X, Y = np.meshgrid( np.linspace(-np.pi,np.pi,51),\
np.linspace(-np.pi,np.pi,51) )
Z = np.sin( X*Y )
fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1,
projection='3d'
)
ax1.plot_surface(X,Y,Z)
ax2 = fig.add_subplot(1, 2, 2)
CS = ax2.contour(X, Y, Z, np.linspace(-1.,1.,10))
ax2.clabel( CS, fontsize=9, inline=1)
plt.show()
P
Y
L
A
B
/
M
A
T
P
L
O
T
L
I
B
It is possible to have an interface (IDE) similar to that of Matlab through Spyder
P
Y
L
A
B
/
M
A
T
P
L
O
T
L
I
B
S
C
I
P
Y
S
C
I
P
Y
Modules of SciPy (reference here:
https://docs.scipy.org/doc/scipy/reference/
)
A
S
S
I
G
N
M
E
N
T
#
1
Create a class to analyze a single degree of freedom system.
The class should initialize an instance with values of mass, stiffness and damping.
In the constructor method, the natural period and pulsation* should be evaluated
and stored in an attribute.
There must be a method to obtain the response under arbitrary loading given as an
array of load values at constant time step; the method should allow to choice
between “central difference” and “Newmark” integration method.
The results of each integration must be saved to a dictionary with an assigned
name, and a method to plot the responses to various loads/integrations should be
present.
Implement the __repr__ method and the comparison method (in terms of natural
period).
*
The International Doctoral Program in Civil and Environmental Engineering at Università degli Studi di Perugia offers a course focused on Python for numerical computing and the development of scientific applications. The program covers interactive command-line interfaces, executing scripts from the OS shell, passing arguments, and using the SciPy ecosystem. It introduces the use of NumPy arrays for multi-dimensional data manipulation, element redistribution among dimensions, and common operations like slicing and iteration. The course provides a practical approach to leveraging Python for advanced engineering applications.
Download Presentation
Please find below an Image/Link to download the presentation.
The content on the website is provided AS IS for your information and personal use only. It may not be sold, licensed, or shared on other websites without obtaining consent from the author. Download presentation by click this link. If you encounter any issues during the download, it is possible that the publisher has removed the file from their server.
E N D
Presentation Transcript
UNIVERSIT DEGLI STUDI DI PERUGIA International Doctoral Program in Civil And Environmental Engineering PYTHON FOR NUMERICAL COMPUTING AND DEVELOPMENT OF SCIENTIFIC APPLICATION Prof. Federico Cluni A. Y. 2020/21 Lesson #3 September, 21th 2021
BASE INTERACTIVITY Command line interface interactivity It is possible to build a script which interact with the user asking for specific data def radici(a,b,c): delta = b**2-4*a*c x1 = (-b-delta**0.5)/(2*a) x2 = (-b+delta**0.5)/(2*a) return x1, x2 a = float(input('Inserisci a: ')) b = float(input('Inserisci b: ')) c = float(input('Inserisci c: ')) x1, x2 = radici(a,b,c) print('x1 = '+str(x1)) print('x2 = '+str(x2)) N.B input read a string from the command line
BASE INTERACTIVITY Command line interface interactivity Moreover, it is also possible to execute a script from the OS shell, and eventually to pass arguments directly from the shell radici.py import sys def radici(a,b,c): delta = b**2-4*a*c x1 = (-b-delta**0.5)/(2*a) x2 = (-b+delta**0.5)/(2*a) return x1, x2 if __name__=='__main__': if len(sys.argv)>1: a = float(sys.argv[1]) b = float(sys.argv[2]) c = float(sys.argv[3]) else: a = float(input('Inserisci a: ')) b = float(input('Inserisci b: ')) c = float(input('Inserisci c: ')) x1, x2 = radici(a,b,c) print('x1 = '+str(x1)) print('x2 = '+str(x2)) > python radici.py > python radici.py 1 1 -2
NUMPY ARRAY < type ndarray > It is a multi-dimensional (may have more than one index) with homogeneous elements and fixed dimension. There are constructors for some typical cases: import numpy as np a = np.zeros(4) a = np.ones((2,3,4)) a = np.arange(-5.,5.,1.) a = np.linspace(-5.,5.,100) N.B this is the standard way of importing numpy Array can also be generated through list and tuples: l = [1., 2., 3. ] a = np.array(l) Remark As usual with object, b=a create an alias and not a copy. To obtain a copy: b = a.copy()
NUMPY ARRAY < type ndarray > It is possible to redistribute the elements among dimension. import numpy.random as rd a = rd.randint(-10,10,(2,3,4)) a.shape = (2,12) N.B numpy.random module is more powerful than random or a.reshape (2,12) Usual slicing operation to extract portions of the array: a[indexes dim. 0, indexes dim. 1, :] Note that slicing produce views of the array, i.e. a[:,:] does not produce a new array! To take all elements of one dimension, use : To access the element at place i a.item(i) The 1-dimensional representation (as a view) is given by a.ravel()
NUMPY ARRAY < type ndarray > To iterate on the dimension for i in range(a.shape[0]): for j in range(a.shape[1]): ... The usual operations +, -, /, *, ** are possible, in such a case the operation is applied to all elements in the array b = 2*a Some useful methods a.max() a.min() a.mean() a.std() a.argmax() a.argmin()
NUMPY MATRIX < type matrix > To define matrices: m = matrix( list/array ) or: m = mat( list/array ) Usual matrix operations are defined (in particular, * is matrix multiplication). To perform matrix multiplication using arrays, use @. In old Python versions, to matrix multiply two array a and b: c=dot(a,b).
NUMPY The Rosetta stone for Matlab users https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
NUMPY - VECTORIALIZATION An array can be used as the argument of a function. Anyway, if there is a if condition, the condition can not be evaluated: def funz(x): if x < 0: return -1 else: return 2*x In such cases, the command where can be used: def funz(x): return np.where( x < 0, -1, 2*x) It is worth noting that this last function can be called with a scalar as an argument. Otherwise funzv = np.vectorize(funz)
NUMPY - USEFUL SUBMODULES/FUNCTIONS There are several submodules in Numpy, for example random, fft, linalg To import a submodule import numpy.random as rd import numpy.linalg as la A = rd.randint(-10,10,(3,3)) b = rd.randint(-10,10,3) if not np.isclose(la.det(A),0,atol=1e-12): c = la.solve(A,b) c1 = np.dot(la.inv(A),b) w, v = la.eig( A ) # solution of A@v = wv
NUMPY - USEFUL SUBMODULES/FUNCTIONS To write/read to a file (in text format) import numpy as np x = np.random.uniform(size=(10,3)) np.savetxt('dati.txt',x,'%7.4f',';') y = np.loadtxt('dati.txt ,delimiter=';')
PYLAB/MATPLOTLIB To create graphs MatLab-like. import numpy as np import matplotlib import matplotlib.pyplot as plt x = np.linspace(0.,2*np.pi, 101) y = np.sin(x) ys = x-x**3/6+x**5/120-x**7/5040+x**9/362880 plt.plot( x, y, 'bo-', label='true') plt.plot( x, ys, 'r-', label='serie 5 terms') plt.xlabel('$t$') plt.ylabel('$\sin(t^1)$') plt.title('Graph of $\sin(t)$') plt.legend() plt.show() To see several examples of use https://matplotlib.org/gallery/index.html
PYLAB/MATPLOTLIB To plot 3d graphs import numpy as np import matplotlib import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D X, Y = np.meshgrid( np.linspace(-np.pi,np.pi,51),\ np.linspace(-np.pi,np.pi,51) ) Z = np.sin( X*Y ) fig = plt.figure() ax1 = fig.add_subplot(1, 2, 1, projection='3d') ax1.plot_surface(X,Y,Z) ax2 = fig.add_subplot(1, 2, 2) CS = ax2.contour(X, Y, Z, np.linspace(-1.,1.,10)) ax2.clabel( CS, fontsize=9, inline=1) plt.show()
PYLAB/MATPLOTLIB It is possible to have an interface (IDE) similar to that of Matlab through Spyder
SCIPY Modules of SciPy (reference here: https://docs.scipy.org/doc/scipy/reference/) cluster Clustering algorithms constants Physical and mathematical constants fftpack Fast Fourier Transform routines integrate Integration and ordinary differential equation solvers interpolate Interpolation and smoothing splines io Input and Output linalg Linear algebra ndimage N-dimensional image processing odr Orthogonal distance regression optimize Optimization and root-finding routines signal Signal processing sparse Sparse matrices and associated routines spatial Spatial data structures and algorithms special Special functions stats Statistical distributions and functions
ASSIGNMENT #1 Create a class to analyze a single degree of freedom system. The class should initialize an instance with values of mass, stiffness and damping. In the constructor method, the natural period and pulsation* should be evaluated and stored in an attribute. There must be a method to obtain the response under arbitrary loading given as an array of load values at constant time step; the method should allow to choice between central difference and Newmark integration method. The results of each integration must be saved to a dictionary with an assigned name, and a method to plot the responses to various loads/integrations should be present. Implement the __repr__ method and the comparison method (in terms of natural period). * 2 T m = = , T 1 1 k 1
