ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 17.04.2021
Просмотров: 3749
Скачиваний: 3
Introduction to Python for Science, Release 0.9.23
45
ax4
.
axhline(color
=
"grey"
, ls
=
"--"
, zorder
=-
1
)
46
ax4
.
set_xlim(
-
15
,
4
)
47
ax4
.
set_ylim(
-
0.5
,
0.6
)
48
ax4
.
text(
0.5
,
0.95
,
’Airy’
, ha
=
’center’
, va
=
’top’
,
49
transform
=
ax4
.
transAxes)
50
51
# Legendre polynomials
52
x
=
np
.
linspace(
-
1
,
1
,
256
)
53
lp0
=
np
.
polyval(sp
.
special
.
legendre(
0
),x)
54
lp1
=
np
.
polyval(sp
.
special
.
legendre(
1
),x)
55
lp2
=
np
.
polyval(sp
.
special
.
legendre(
2
),x)
56
lp3
=
np
.
polyval(sp
.
special
.
legendre(
3
),x)
57
ax5
=
fig
.
add_subplot(
325
)
58
ax5
.
plot(x,lp0, x,lp1, x,lp2, x,lp3)
59
ax5
.
axhline(color
=
"grey"
, ls
=
"--"
, zorder
=-
1
)
60
ax5
.
set_ylim(
-
1
,
1.1
)
61
ax5
.
text(
0.5
,
0.9
,
’Legendre’
, ha
=
’center’
, va
=
’top’
,
62
transform
=
ax5
.
transAxes)
63
64
# Laguerre polynomials
65
x
=
np
.
linspace(
-
5
,
8
,
256
)
66
lg0
=
np
.
polyval(sp
.
special
.
laguerre(
0
),x)
67
lg1
=
np
.
polyval(sp
.
special
.
laguerre(
1
),x)
68
lg2
=
np
.
polyval(sp
.
special
.
laguerre(
2
),x)
69
lg3
=
np
.
polyval(sp
.
special
.
laguerre(
3
),x)
70
ax6
=
fig
.
add_subplot(
326
)
71
ax6
.
plot(x,lg0, x,lg1, x,lg2, x,lg3)
72
ax6
.
axhline(color
=
"grey"
, ls
=
"--"
, zorder
=-
1
)
73
ax6
.
axvline(color
=
"grey"
, ls
=
"--"
, zorder
=-
1
)
74
ax6
.
set_xlim(
-
5
,
8
)
75
ax6
.
set_ylim(
-
5
,
10
)
76
ax6
.
text(
0.5
,
0.9
,
’Laguerre’
, ha
=
’center’
, va
=
’top’
,
77
transform
=
ax6
.
transAxes)
78
79
plt
.
show()
The arguments of the different functions depend, of course, on the nature of the particular
function. For example, the first argument of the two types of Bessel functions called
in lines 10-13 is the so-called
order
of the Bessel function, and the second argument is
the independent variable. The Gamma and Error functions take one argument each and
produce one output. The Airy function takes only one input argument, but returns four
outputs, which correspond the two Airy functions, normally designated
Ai(
x
)
and
Bi(
x
)
,
and their derivatives
Ai
0
(
x
)
and
Bi
0
(
x
)
. The plot shows only
Ai(
x
)
and
Bi(
x
)
.
160
Chapter 9. Numerical Routines: SciPy and NumPy
Introduction to Python for Science, Release 0.9.23
The polynomial functions shown have a special syntax that uses NumPy’s
polyval
function for generating polynomials. If
p
is a list or array of
N
numbers and
x
is an array,
then
polyval(p, x)
=
p[
0
]
*
x
**
(N
-
1
)
+
p[
1
]
*
x
**
(N
-
2
)
+ ... +
p[N
-
2
]
*
x
+
p[N
-
1
]
For example, if
p = [2.0, 5.0, 1.0]
,
polyval
generates the following
quadratic polynomial:
polyval(p, x)
=
2.0
*
x
**
2
+
5.0
*
x
+
1.0
SciPy’s
special.legendre(n)
and
special.laguerre(n)
functions output
the coefficients
p
needed in
polyval
to produce the
n
th
-order Legendre and Laguerre
polynomials, respectively. The
special
library of SciPy has functions that specify
many other polynomial functions in this same way.
9.2 Linear algebra
Python’s mathematical libraries, NumPy and SciPy, have extensive tools for numerically
solving problems in linear algebra. Here we focus on two problems that arise commonly
in scientific and engineering settings: (1) solving a system of linear equations and (2)
eigenvalue problems. In addition, we also show how to perform a number of other basic
computations, such as finding the determinant of a matrix, matrix inversion, and
LU
decomposition. The SciPy package for linear algebra is called
scipy.linalg
.
9.2.1 Basic computations in linear algebra
SciPy has a number of routines for performing basic operations with matrices. The deter-
minant of a matrix is computed using the
scipy.linalg.det
function:
In [1]:
import
scipy.linalg
In [2]:
a
=
array([[
-
2
,
3
], [
4
,
5
]])
In [3]:
a
Out[4]:
array([[
-
2
,
3
],
[ 4,
5]])
In [5]:
scipy
.
linalg
.
det(a)
Out[5]:
-
22.0
The inverse of a matrix is computed using the
scipy.linalg.inv
function, while the
product of two matrices is calculated using the NumPy
dot
function:
9.2. Linear algebra
161
Introduction to Python for Science, Release 0.9.23
In [6]:
b
=
scipy
.
linalg
.
inv(a)
In [6]:
b
Out[6]:
array([[
-
0.22727273
,
0.13636364
],
[ 0.18181818,
0.09090909]])
In [7]:
dot(a,b)
Out[7]:
array([[
1.
,
0.
],
[ 0.,
1.]])
9.2.2 Solving systems of linear equations
Solving systems of equations is nearly as simple as constructing a coefficient matrix and
a column vector. Suppose you have the following system of linear equations to solve:
2
x
1
+ 4
x
2
+ 6
x
3
= 4
x
1
−
3
x
2
−
9
x
3
=
−
11
8
x
1
+ 5
x
2
−
7
x
3
= 1
The first task is to recast this set of equations as a matrix equation of the form
A
x
=
b
.
In this case, we have:
A
=
2
4
6
1
−
3
−
9
8
5
−
7
,
x
=
x
1
x
2
x
3
,
b
=
4
−
11
1
.
Next we construct the array
A
and vector
b
as NumPy arrays:
In [8]:
A
=
array([[
2
,
4
,
6
], [
1
,
-
3
,
-
9
], [
8
,
5
,
-
7
]])
In [9]:
b
=
array([
4
,
-
11
,
2
])
Finally we use the SciPy function
scipy.linalg.solve
to find
x
1
,
x
2
, and
x
3
.
In [10]:
scipy
.
linalg
.
solve(A,b)
Out[10]:
array([
-
8.91304348
,
10.2173913
,
-
3.17391304
])
which gives the results:
x
1
=
−
8
.
91304348
,
x
2
= 10
.
2173913
, and
x
3
=
−
3
.
17391304
.
Of course, you can get the same answer by noting that
x
=
A
−
1
b
. Following this ap-
proach, we can use the
scipy.linalg.inv
introduced in the previous section:
Ainv = scipy.linalg.inv(A)
In [10]:
dot(Ainv, b)
Out[10]:
array([
-
8.91304348
,
10.2173913
,
-
3.17391304
])
162
Chapter 9. Numerical Routines: SciPy and NumPy
Introduction to Python for Science, Release 0.9.23
which is the same answer we obtained using
scipy.linalg.solve
.
Using
scipy.linalg.solve
is numerically more stable and a faster than using
x
=
A
−
1
b
,
so it is the preferred method for solving systems of equations.
You might wonder what happens if the system of equations are not all linearly indepen-
dent. For example if the matrix
A
is given by
A
=
2
4
6
1
−
3
−
9
1
2
3
where the third row is a multiple of the first row. Let’s try it out and see what happens.
First we change the bottom row of the matrix
A
and then try to solve the system as we did
before.
In [11]:
A[
2
]
=
array([
1
,
2
,
3
])
In [12]:
A
Out[12]:
array([[
2
,
4
,
6
],
[ 1, -3, -9],
[ 1,
2,
3]])
In [13]:
scipy
.
linalg
.
solve(A,b)
LinAlgError: Singular matrix
In [14]:
Ainv
=
scipy
.
linalg
.
inv(A)
LinAlgError: Singular matrix
Whether we use
scipy.linalg.solve
or
scipy.linalg.inv
, SciPy raises an
error because the matrix is singular.
9.2.3 Eigenvalue problems
One of the most common problems in science and engineering is the eigenvalue problem,
which in matrix form is written as
A
x
=
λ
x
where
A
is a square matrix,
x
is a column vector, and
λ
is a scalar (number). Given
the matrix
A
, the problem is to find the set of eigenvectors
x
and their corresponding
eigenvalues
λ
that solve this equation.
We can solve eigenvalue equations like this using
scipy.linalg.eig
. the outputs of
this function is an array whose entries are the eigenvalues and a matrix whose rows are the
9.2. Linear algebra
163
Introduction to Python for Science, Release 0.9.23
eigenvectors. Let’s return to the matrix we were using previously and find its eigenvalues
and eigenvectors.
A = array([[2, 4, 6],[1, -3, -9],[8, 5, -7]])
In [15]:
A
Out[15]:
array([[
2
,
4
,
6
],
[ 1, -3, -9],
[ 8,
5, -7]])
In [16]:
lam, evec
=
scipy
.
linalg
.
eig(A)
In [17]:
lam
Out[17]:
array([
2.40995356
+
0.j
,
-
8.03416016
+
0.j
,
-2.37579340+0.j])
In [18]:
evec
Out[18]:
array([[
-
0.77167559
,
-
0.52633654
,
0.57513303
],
[ 0.50360249,
0.76565448, -0.80920669],
[-0.38846018,
0.36978786,
0.12002724]])
The first eigenvalue and its corresponding eigenvector are given by
In [19]:
lam[
0
]
Out[19]:
(
2.4099535647625494
+
0j
)
In [20]:
evec[:,
0
]
Out[20]:
array([
-
0.77167559
,
0.50360249
,
-
0.38846018
])
We can check that they satisfy the
A
x
=
λ
x
:
In [21]:
dot(A,evec[:,
0
])
Out[21]:
array([
-
1.85970234
,
1.21365861
,
-
0.93617101
])
In [22]:
lam[
0
]
*
evec[:,
0
]
Out[22]:
array([
-
1.85970234
+
0.j
,
1.21365861
+
0.j
,
-0.93617101+0.j])
Thus we see by direct substitution that the left and right sides of
A
x
=
λ
x
:
are equal. In
general, the eigenvalues can be complex, so their values are reported as complex numbers.
164
Chapter 9. Numerical Routines: SciPy and NumPy