Python科学计算 第3版pdf python科学计算第三版

本系列是使用Ubuntu 64位操作系统上，anaconda环境下，采用jupyter notebook 实现python3科学计算的文章，应用于青岛科技大学信息学院科学计算与数据挖掘等多个课程中python教学。需要安装的模块，numpy ，matplotlib，scipy，sympy，pandas等科学计算常用模块。由于匆忙成稿，个人能力所限，显得粗糙，并且有几个地方是错误的。主要参考了Numerical Python 2nd Edition

科学计算包括数值计算和符号计算两种计算 。计算机能够对数值进行一系列运算是人所共知的事，但计算机也能够对含未知量的式子直接进行推导、演算则并不是人人皆知。数值计算和符号计算本来应该是并存的两种计算，是计算的平行的两个部分，决不能厚此薄彼，因此这两种计算都是一样重要的。利用计算机对一个函数进行求导、积分，这早己成为事实。

import sympy

sympy.init_printing()#调用及时显示系统。但对于jupyter似乎不需要。

x = sympy.Symbol("x")#创建符号

符号名是一个字符串，可以选择包含类似Latex标记

使符号名在ipython的显示系统中显示良好。

符号对象的名称在创建时设置。符号可以采用不同方式创建，

例如，使用sympy.symbol、sympy.symbols和sympy。

通常需要将sympy符号与python变量关联起来

同一个名称或与符号名称紧密对应的名称。例如,

要创建一个名为x的符号，并将其绑定到同名的python变量，

我们可以使用symbol类的构造函数并传递包含该符号的字符串

名称作为第一个参数

y = sympy.Symbol("y", real=True)

y.is_real

True

x.is_real is None

True

另一方面，如果我们要使用is real来查询先前定义的符号

x，它没有显式地指定为real，因此可以同时表示real和

非实变量，结果是没有.但如果被显式定义，则为false

y.is_real is None

False

sympy.Symbol("z", imaginary=True).is_real

False

x = sympy.Symbol("x")

y = sympy.Symbol("y", positive=True)

sympy.sqrt(x ** 2)

$\displaystyle \sqrt{x^{2}}$

sympy.sqrt(y ** 2)

$\displaystyle y$

#以上运算可以清晰的发现定义的清晰度，可以影响运算结果。

n1 = sympy.Symbol("n")

n2 = sympy.Symbol("n", integer=True)

n3 = sympy.Symbol("n", odd=True)

import numpy as np

sympy.cos(n1 * np.pi)

$\displaystyle \cos{\left(3.14159265358979 n \right)}$

sympy.cos(n2 * np.pi)

$\displaystyle \cos{\left(3.14159265358979 n \right)}$

sympy.cos(n3 *np.pi)

$\displaystyle \cos{\left(3.14159265358979 n \right)}$

from sympy import I, pi, oo #然而在sympy中不建议采用sympy之外的常数

sympy.cos(n1 * pi)

$\displaystyle \cos{\left(\pi n \right)}$

sympy.cos(n2 * pi)

$\displaystyle \left(-1\right)^{n}$

sympy.cos(n3 *pi)

$\displaystyle -1$

可以看到采用sympy内部的常数，虚数，π，无穷大，更利于公式简化推导。

oo

$\displaystyle \infty$

pi

$\displaystyle \pi$

sympy.E

$\displaystyle e$

sympy.EulerGamma

$\displaystyle \gamma$

sympy.I

$\displaystyle i$

sympy.oo

$\displaystyle \infty$

创建符号时，符号的属性可以用下列属性标定：

real ,imaginary 属性判定 is_real, is_imaginary

positive, negative 属性判定 is_positive, is_negative

integer 属性判定 is_integer

odd, even is_odd, 属性判定 is_even

prime 属性判定 is_prime

finite, infinite 属性判定 is_finite, is_infinite

a, b, c = sympy.symbols("a, b, c", negative=True)

d, e, f = sympy.symbols("d, e, f", positive=True)

i = sympy.Integer(19)

type(i)

sympy.core.numbers.Integer

i.is_Integer, i.is_real, i.is_odd

(True, True, True)

"%.25f" % 0.3  # create a string representation with 25 decimals，把0.3表示成有25个小数的浮点数。

'0.2999999999999999888977698'

sympy.Float(0.3, 25)

$\displaystyle 0.2999999999999999888977698$

sympy.Float('0.3', 25)#请同学们注意的是，这个才是精确的表示浮点数0.3

$\displaystyle 0.3$

但是，请注意，要正确地将0.3表示为浮点对象，必须

从字符串“0.3”而不是python float 0.3初始化它，后者已经包含

浮点错误。

r1 = sympy.Rational(2, 3)#分式的表示方式

r1

$\displaystyle \frac{2}{3}$

5/6

$\displaystyle 0.8333333333333334$

x/y

$\displaystyle \frac{x}{y}$

定义符号函数

x,y,z = sympy.symbols("x,y,z")

f = sympy.Function("f")

type(f)

sympy.core.function.UndefinedFunction

f(x)

$\displaystyle f{\left(x \right)}$

g = sympy.Function("g")(x,y,z)

g

$\displaystyle g{\left(x,y,z \right)}$

g.free_symbols

$\displaystyle \left\{x, y, z\right\}$

g

$\displaystyle g{\left(x,y,z \right)}$

定义函数符号，下面可以定义函数的表达式

h = sympy.Lambda(x,x**2+2*x)

h

$\displaystyle \left( x \mapsto x^{2} + 2 x \right)$

#lambda函数定义一类函数只有当你的情况完全满足这四个标准时, 我才会说你可以使用 lambda 表达式:
所谓的匿名函数

你所要做的操作是不重要的:函数不值得一个名称

使用 lambda 表达式比你所能想到的函数名称让代码更容易理解

你很确定还没有一个函数能满足你的需求

你团队的每个人都了解 lambda 表达式, 并且都同意使用它们

如果上述四条中的任何一条都不符合你的情况, 我建议用 def 来写一个新的函数, (如果可能)接受一个在 Python 中已经存在且能满足你需求的函数.

def f(x):
    return x**2+2*x

f(x)

$\displaystyle x^{2} + 2 x$

#很多时候，我们并不需要函数表示，而是需要如下的表达式的(柄,感觉叫句柄有点意思)
expre = 1+x**2+5*x**3

expre

$\displaystyle 5 x^{3} + x^{2} + 1$

#expre,有自己的完善的参数体系，这个体系呈现树形结构。

expre.args

$\displaystyle \left( 1, \ x^{2}, \ 5 x^{3}\right)$

expre.args[2]

$\displaystyle 5 x^{3}$

expre.args[2].args[1]

$\displaystyle x^{3}$

expre.args[2].args[0]

$\displaystyle 5$

#化简操作，是所有表达式运算过程中最常用的函数

expr = 2*(x**2-x) - x*(x+1)

expr

$\displaystyle 2 x^{2} - x \left(x + 1\right) - 2 x$

sympy.simplify(expr)

$\displaystyle x \left(x - 3\right)$

expr2 =2*sympy.cos(x)*sympy.sin(x)

expr2

$\displaystyle 2 \sin{\left(x \right)} \cos{\left(x \right)}$

sympy.simplify(expr2)

$\displaystyle \sin{\left(2 x \right)}$

化简成功的关键是对所有的函数和常数采用sympy定义的函数或常数。

化简函数

sympy.simplify 调用不同的方法得到化简结果

sympy.trigsimp 尝试三角恒等式化简结果

sympy.powsimp 用幂指数化简

sympy.compsimp 通过组合化简多项式

sympy.ratsimp 提取公分母

下面是展开函数的使用说明：

expr=(x+4)*(x+5)

sympy.expand(expr)

$\displaystyle x^{2} + 9 x + 20$

sympy.sin(x+y).expand()

$\displaystyle \sin{\left(x + y \right)}$

sympy.sin(x+y).expand(trig= True)

$\displaystyle \sin{\left(x \right)} \cos{\left(y \right)} + \sin{\left(y \right)} \cos{\left(x \right)}$

a, b = sympy.symbols("a, b", positive=True)# 使用化简必须保证，所有的数学符号，函数都是sympy产生的

sympy.log(a * b).expand(log=True)

$\displaystyle \log{\left(a \right)} + \log{\left(b \right)}$

sympy.exp(I*a + b).expand(complex=True)

$\displaystyle i e^{b} \sin{\left(a \right)} + e^{b} \cos{\left(a \right)}$

sympy.expand((a * b)**x, power_base=True)

$\displaystyle a^{x} b^{x}$

sympy.exp((a-b)*x).expand(power_exp=True)

$\displaystyle e^{a x} e^{- b x}$

因子、汇集和组合

sympy.factor(x**2 - 1)

$\displaystyle \left(x - 1\right) \left(x + 1\right)$

sympy.factor(x * sympy.cos(y) + sympy.sin(z) * x)

$\displaystyle x \left(\sin{\left(z \right)} + \cos{\left(y \right)}\right)$

sympy.logcombine(sympy.log(a) - sympy.log(b))

$\displaystyle \log{\left(\frac{a}{b} \right)}$

sympy.trigsimp, sympy.powsimp, and sympy.logcombine,都可以产生简单的化简

expr = x + y + x**2 * y * z

expr.collect(x) #按照x的指数大小排列

$\displaystyle x^{2} y z + x + y$

expr = sympy.cos(x + y) + sympy.sin(x - y)

expr.expand(trig=True).collect([sympy.cos(x), 
                                 sympy.sin(x)]).collect(sympy.cos(y) - sympy.sin(y))

$\displaystyle \left(\sin{\left(x \right)} + \cos{\left(x \right)}\right) \left(- \sin{\left(y \right)} + \cos{\left(y \right)}\right)$

expr.expand(trig=True).
collect([sympy.cos(x), sympy.sin(x)]).collect(sympy.cos(y) - sympy.sin(y))

File "<ipython-input-96-fefb96c05190>", line 1
    expr.expand(trig=True).
                           ^
SyntaxError: invalid syntax

expr.expand(trig=True).collect([
    sympy.cos(x), sympy.sin(x)]).collect(sympy.cos(y) - sympy.sin(y))

$\displaystyle \left(\sin{\left(x \right)} + \cos{\left(x \right)}\right) \left(- \sin{\left(y \right)} + \cos{\left(y \right)}\right)$

事实上，在那个地方回车是有讲究的，如果愿意全部试试。

sympy.apart(1/(x**2 + 3*x + 2), x)#非常有用的分解

$\displaystyle - \frac{1}{x + 2} + \frac{1}{x + 1}$

sympy.together(1 / (y * x + y) + 1 / (1+x))

$\displaystyle \frac{y + 1}{y \left(x + 1\right)}$

sympy.cancel(y / (y * x + y)) #约分化简

$\displaystyle \frac{1}{x + 1}$

sympy.cancel(y / (y * x + 4))

$\displaystyle \frac{y}{x y + 4}$

sympy.cancel(y / (y * x + 4*y**2))

$\displaystyle \frac{1}{x + 4 y}$

代换 Substitutions

(x + y).subs(x, y)

$\displaystyle 2 y$

sympy.sin(x * sympy.exp(x)).subs(x, y)

$\displaystyle \sin{\left(y e^{y} \right)}$

sympy.sin(x * z).subs({z: sympy.exp(y), x: y, sympy.sin: sympy.cos})

$\displaystyle \cos{\left(y e^{y} \right)}$

当然也可以带入值

expr = x * y + z**2 *x

values = {x: 1.25, y: 0.4, z: 3.2}

expr.subs(values)

$\displaystyle 13.3$

赋值过程中，追求精度，必须要采用，sympy.N()函数，或者用evalf()函数、

sympy.N(pi, 50)

$\displaystyle 3.1415926535897932384626433832795028841971693993751$

(x + 1/pi).evalf(10)

$\displaystyle x + 0.3183098862$

f = sympy.Function('f')(x)

sympy.diff(f, x)               # equivalent to f.diff(x)

$\displaystyle \frac{d}{d x} f{\left(x \right)}$

sympy.diff(f, x, x)

$\displaystyle \frac{d^{2}}{d x^{2}} f{\left(x \right)}$

sympy.diff(f, x, 3)   # equivalent to sympy.diff(f, x, x, x)

$\displaystyle \frac{d^{3}}{d x^{3}} f{\left(x \right)}$

g = sympy.Function('g')(x, y)

g.diff(x, y)          # equivalent to sympy.diff(g, x, y)

$\displaystyle \frac{\partial^{2}}{\partial y\partial x} g{\left(x,y \right)}$

g.diff(x, 3, y, 2)    # equivalent to sympy.diff(g, x, x, x, y, y)

$\displaystyle \frac{\partial^{5}}{\partial y^{2}\partial x^{3}} g{\left(x,y \right)}$

expr = x**4 + x**3 + x**2 + x + 1

expr.diff(x)

$\displaystyle 4 x^{3} + 3 x^{2} + 2 x + 1$

expr = sympy.sin(x * y) * sympy.cos(x / 2)

expr.diff(x)

$\displaystyle y \cos{\left(\frac{x}{2} \right)} \cos{\left(x y \right)} - \frac{\sin{\left(\frac{x}{2} \right)} \sin{\left(x y \right)}}{2}$

expr2 = sympy.special.polynomials.hermite(x, 0)

expr.diff(x).doit()

$\displaystyle \frac{2^{x} \sqrt{\pi} \operatorname{polygamma}{\left(0,\frac{1}{2} - \frac{x}{2} \right)}}{2 \Gamma\left(\frac{1}{2} - \frac{x}{2}\right)} + \frac{2^{x} \sqrt{\pi} \log{\left(2 \right)}}{\Gamma\left(\frac{1}{2} - \frac{x}{2}\right)}$

d = sympy.Derivative(sympy.exp(sympy.cos(x)), x)

d

$\displaystyle \frac{d}{d x} e^{\cos{\left(x \right)}}$

d.doit()#执行上一步运算

$\displaystyle - e^{\cos{\left(x \right)}} \sin{\left(x \right)}$

a, b, x, y = sympy.symbols("a, b, x, y")

f = sympy.Function("f")(x)

sympy.integrate(f)

$\displaystyle \int f{\left(x \right)}\, dx$

sympy.integrate(f, (x, a, b))

$\displaystyle \int\limits_{a}^{b} f{\left(x \right)}\, dx$

sympy.integrate(sympy.exp(-x**2), (x, 0, oo))

$\displaystyle \frac{\sqrt{\pi}}{2}$

级数

x, y = sympy.symbols("x, y")

f = sympy.Function("f")(x)

sympy.series(f, x)

KaTeX parse error: Undefined control sequence: \substack at position 87: …ght)} \right|_{\̲s̲u̲b̲s̲t̲a̲c̲k̲{ x=0 }} + \fra…

x0 = sympy.Symbol("{x_0}")

f.series(x, x0, n=3)

KaTeX parse error: Undefined control sequence: \substack at position 124: …ght)} \right|_{\̲s̲u̲b̲s̲t̲a̲c̲k̲{ \xi_{1}={x_0}…

f.series(x, x0, n=2).removeO() #移除掉截断项，参与后续运算。

KaTeX parse error: Undefined control sequence: \substack at position 99: …ght)} \right|_{\̲s̲u̲b̲s̲t̲a̲c̲k̲{ \xi_{1}={x_0}…

f.series(x, x0, n=3).removeO()

KaTeX parse error: Undefined control sequence: \substack at position 117: …ght)} \right|_{\̲s̲u̲b̲s̲t̲a̲c̲k̲{ \xi_{1}={x_0}…

sympy.cos(x).series()

$\displaystyle 1 - \frac{x^{2}}{2} + \frac{x^{4}}{24} + O\left(x^{6}\right)$

(1/(1+x)).series()

$\displaystyle 1 - x + x^{2} - x^{3} + x^{4} - x^{5} + O\left(x^{6}\right)$

expr = sympy.cos(x) / (1 + sympy.sin(x * y))

expr.series(x, n=4)

$\displaystyle 1 - x y + x^{2} \left(y^{2} - \frac{1}{2}\right) + x^{3} \left(- \frac{5 y^{3}}{6} + \frac{y}{2}\right) + O\left(x^{4}\right)$

expr.series(y, n=4)

$\displaystyle \cos{\left(x \right)} - x y \cos{\left(x \right)} + x^{2} y^{2} \cos{\left(x \right)} - \frac{5 x^{3} y^{3} \cos{\left(x \right)}}{6} + O\left(y^{4}\right)$

sympy.limit(sympy.sin(x) / x, x, 0) #极限

$\displaystyle 1$

f = sympy.Function('f')

x, h = sympy.symbols("x, h")

diff_limit = (f(x + h) - f(x))/h

sympy.limit(diff_limit.subs(f, sympy.cos), h, 0)

$\displaystyle - \sin{\left(x \right)}$

sympy.limit(diff_limit.subs(f, sympy.sin), h, 0)

$\displaystyle \cos{\left(x \right)}$

以上代码用于说名sin函数的极限由来，没有太强的意义

expr = (x**2 - 3*x) / (2*x - 2)

p = sympy.limit(expr/x, x, sympy.oo)

q = sympy.limit(expr - p*x, x, sympy.oo)

p,q

$\displaystyle \left( \frac{1}{2}, \ -1\right)$

n = sympy.symbols("n", integer=True)

x = sympy.Sum(1/(n**2), (n, 1, oo))

n

$\displaystyle n$

x

$\displaystyle \sum_{n=1}^{\infty} \frac{1}{n^{2}}$

x.doit() #有的函数必须用doit（）函数来完成，最终的运算

$\displaystyle \frac{\pi^{2}}{6}$

x = sympy.Symbol("x")

sympy.Sum((x)**n/(sympy.factorial(n)), (n, 1, oo)).doit().simplify()

$\displaystyle e^{x} - 1$

wz = sympy.Sum((x)**n/(sympy.factorial(n)), (n, 1, oo))

wz.doit()

$\displaystyle x \left(\frac{e^{x}}{x} - \frac{1}{x}\right)$

wz.doit().simplify()

$\displaystyle e^{x} - 1$

涉及到各类化简，再加上doit（）命令，在函数化简的路上，要多尝试！

方程

x = sympy.Symbol("x")

sympy.solve(x**2 + 2*x - 3)

$\displaystyle \left[ -3, \ 1\right]$

a, b, c = sympy.symbols("a, b, c")

sympy.solve(a * x**2 + b * x + c, x)

$\displaystyle \left[ \frac{- b + \sqrt{- 4 a c + b^{2}}}{2 a}, \ - \frac{b + \sqrt{- 4 a c + b^{2}}}{2 a}\right]$

eq1 = x + 2 * y - 1

eq2 = x - y + 1

sympy.solve([eq1, eq2], [x, y], dict=True)

$\displaystyle \left[ \left\{ x : - \frac{1}{3}, \ y : \frac{2}{3}\right\}\right]$

#也可以获取结果
sols = sympy.solve([eq1, eq2], [x, y], dict=True)

[eq1.subs(sol).simplify() == 0 and eq2.subs(sol).simplify() == 0 for sol in sols]

[True]

线性代数

sympy.Matrix([1,2])

$\displaystyle \left[\begin{matrix}1\\2\end{matrix}\right]$

np.array([1,2])

array([1, 2])

a, b, c, d = sympy.symbols("a, b, c, d")

M = sympy.Matrix([[a, b], [c, d]])

M

$\displaystyle \left[\begin{matrix}a & b\\c & d\end{matrix}\right]$

M * M

$\displaystyle \left[\begin{matrix}a^{2} + b c & a b + b d\\a c + c d & b c + d^{2}\end{matrix}\right]$

M.norm()

$\displaystyle \sqrt{\left|{a}\right|^{2} + \left|{b}\right|^{2} + \left|{c}\right|^{2} + \left|{d}\right|^{2}}$

M.rank()

$\displaystyle 2$

M.trace()

$\displaystyle a + d$

M.transpose()

$\displaystyle \left[\begin{matrix}a & c\\b & d\end{matrix}\right]$

M^T

---------------------------------------------------------------------------

NameError                                 Traceback (most recent call last)

<ipython-input-193-0be5c4114203> in <module>
----> 1 M^T


NameError: name 'T' is not defined

M.T

$\displaystyle \left[\begin{matrix}a & c\\b & d\end{matrix}\right]$

M.H

$\displaystyle \left[\begin{matrix}\overline{a} & \overline{c}\\\overline{b} & \overline{d}\end{matrix}\right]$

transpose/T 转置
adjoint/H 转置加共轭
trace 求迹
det 行列式
inv 矩阵的逆
LUdecomposition LU分解
LUsolve LU求解方程
QRdecomposition 分解
diagonalize 对角化 D = P^(-1) MP
norm 范数
nullspace Compute a set of vectors that span the null space of a Matrix.
rank 计算矩阵的秩
singular_values 奇异值
solve Solve a linear system of equations in the form Mx = b.

example： $ x+py = b_{1} $,$ qx+y = b_{2} $
我们可以把上面的方程用矩阵表示为
$ MX =b $

p, q = sympy.symbols("p, q")

M = sympy.Matrix([[1, p], [q, 1]])

b = sympy.Matrix(sympy.symbols("b_1, b_2"))

x = M.LUsolve(b)

x

$\displaystyle \left[\begin{matrix}b_{1} - \frac{p \left(- b_{1} q + b_{2}\right)}{- p q + 1}\\\frac{- b_{1} q + b_{2}}{- p q + 1}\end{matrix}\right]$

x1 = M.inv() * b

x1

$\displaystyle \left[\begin{matrix}\frac{b_{1}}{- p q + 1} - \frac{b_{2} p}{- p q + 1}\\- \frac{b_{1} q}{- p q + 1} + \frac{b_{2}}{- p q + 1}\end{matrix}\right]$

x ==x1

False

rt = x.T*x1
rt

$\displaystyle \left[\begin{matrix}\frac{\left(- b_{1} + b_{2} p\right) \left(b_{1} - \frac{p \left(- b_{1} q + b_{2}\right)}{- p q + 1}\right)}{p q - 1} + \frac{\left(- b_{1} q + b_{2}\right) \left(b_{1} q - b_{2}\right)}{\left(- p q + 1\right) \left(p q - 1\right)}\end{matrix}\right]$

rt.collect(p)#调用化简函数都会报错，原因在于多变量矩阵没有该类化简的属性

---------------------------------------------------------------------------

AttributeError                            Traceback (most recent call last)

<ipython-input-232-b3d07a6ec5ae> in <module>
----> 1 rt.collect(p)


~/.conda/envs/qutip-env/lib/python3.7/site-packages/sympy/matrices/matrices.py in __getattr__(self, attr)
   2139         else:
   2140             raise AttributeError(
-> 2141                 "%s has no attribute %s." % (self.__class__.__name__, attr))
   2142 
   2143     def __len__(self):


AttributeError: MutableDenseMatrix has no attribute collect.

ics={b[0]:3,b[1]:4,p:7,q:11}
x.subs(ics)

$\displaystyle \left[\begin{matrix}\frac{25}{76}\\\frac{29}{76}\end{matrix}\right]$

b

$\displaystyle \left[\begin{matrix}b_{1}\\b_{2}\end{matrix}\right]$

b[0]

$\displaystyle b_{1}$

b[1]

$\displaystyle b_{2}$

x1.subs(ics)

$\displaystyle \left[\begin{matrix}\frac{25}{76}\\\frac{29}{76}\end{matrix}\right]$

由此可以看到，直接比较是会得到荒谬的答案的，然而带入数据后，发现其实是相等的。

然而，如果解方程mx=b是目标，计算矩阵的逆比执行LU分解的难度更大，

那么使用lu分解更有效。这对于较大的方程系统，尤其重要