Python數學運算與科學運算

常見套件

也可以使用R語言替代Python

安裝

archlinux安裝上面所提到的套件套件

sudo pacman -S python-scipy python-sympy python-numpy python-matplotlib

Numpy使用

匯入numpy

import numpy as np

關於常用的數學函數參見numpy的數學函式表

定義陣列並且運算，numpy的陣列可以相加或相乘內容， 不同於python本身的陣列只是將兩個陣列合而為一。

x1 = np.array([1,2,3,4])
x2 = np.array([1,2,3,4])
x1 + x2
x1 * x2

二維矩陣相乘

x1 = np.array([[2,5],[-1,4]])
x2 = np.array([[0, -2, 0],[1,-3,6]])
x1.dot(x2)

容易混淆的部份是數學上的相乘應該使用array.dot(array)的方式使用， 而不是array1 * array2，下面這個範例可以看出不同。

x1 = np.array([[1,2],[3,4]])
x2 = np.array([[1,2],[3,4]])
x1 * x2
x1.dot(x2)

平均值、變異數、標準差

x = np.array([1, 2, 3, 4])
np.mean(x)
np.var(x)
np.std(x)

求總和，不要和sympy的Sum()搞混

a = np.array([1,2,3,4])
np.sum(a)

Sympy

微分積分

導入

from sympy import *

定義符號

x, y, z = symbols('x y z')

微分與積分

diff(sin(x), x)
integrate(1/x, x)

泰勒級數(任何連續函數都可以用多項式表示)

expr = sin(x)
expr.series(x, 0, 4)

解方程式

x, y = symbols('x y')
e1 = 2 * x + y
e2 = 4 * x + y + 2
solve([e1, e2], [x, y])

代入常數或變數

x, y = symbols('x y')
e1 = 2 * x + y
e1.evalf(subs={x: 1, y:2})
e1.evalf(subs={x: 1, y:x})

注意: evalf()進行計算的行為，代換變數採用subs()

x, y = symbols('x y')
e1 = 2 * x + y
e1.subs({x: 1, y:2})
e1.subs({x: 1, y:x})

左右式

x = symbols('x')
eq = Eq(x * 2 + 1, x ** 2)
eq.lhs
eq.rhs

級數

我們可以透過sympy來求級數，可以使用Sum()這個指令， 切記不要和np.sum()搞混。

可以透過以下方法來求總和，我們求一個1, 2, 3, ..., 10數列總和

x = symbols('x')
f = Sum(x, (x, 0, 10))
f.doit()

在上面我們定義了一個等差級數f，並且透過f.doit()求解， 而(x, 0, 10)的0和10代表開始與結尾。

而sympy也可以求出無窮級數，使用oo作為無限，我們定義一個 e ^ x的泰勒級數，並且採用doit()來產生e ^ x。

x, n = symbols('x n')
f = Sum(x ** n / factorial(n), (n, 0, oo))
f.doit()

執行之後我們可以看出sympy給出了exp(x)

微分方程式求解

先看我們範例程式

x = symbols('x')
y = symbols('y', cls=Function)

diffeq = Eq(y(x).diff(x) + y(x), 5)
ans = dsolve(diffeq, y(x))
print(ans)

第一步，宣告自變數x與應變數y

x = symbols('x')
y = symbols('y', cls=Function)

第二步，定義微分方程式

diffeq = Eq(y(x).diff(x) + y(x), 5)

第三步，求解

ans = dsolve(diffeq, y(x))
print(ans)

進一步解決IVP，第三步請改成以下的程式碼， 表示設定一個初始條件y(0) = 1

ans = dsolve(diffeq, y(x), ics={y(0): 1})
print(ans)

請參見sympy文檔

latex

latex()可以用來產生latex語法，而不必自己編寫排版

Matplitlib

引用Matplitlib與numpy

import matplotlib.pyplot as plt
import numpy as np

畫出圖形，可以利用numpy的linspace()來生成點，區間為0~100並且分成20份

x = np.linspace(0, 100, 20)
y = x ** 2
plt.plot(x, y)

設定圖表參數

plt.xlabel('X')
plt.ylabel('f(x) = y')
plt.title('This is title')

可以加上顯示文字，我們將他放置在20,400的位置

plt.text(20,400,'(20, 400)')

顯示圖表，可能會發現圖形沒有很精準在數值上，可以把linspace()調整精細

plt.show()

另外可以採用savefig()來保存圖示，可以保存為png或svg。

plt.savefig('fig.svg')
plt.savefig('fig.png')

採用clf()清理上次所畫出圖

plt.clf()

顯示網格

plt.grid(True)

圖例

plt.legend(('data1', 'data2', 'data3'))

範例程式

範例程式1

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    y = np.exp(-x) * np.sin(2 * x)
    return y

x = np.linspace(0, 2 * np.pi, 100)
y = f(x)

plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Plot of the Function f(x) = exp(-x) * sin(2x)')
plt.text(4, -0.1, 'Copyright@lin')

plt.show()

範例程式2，繪製兩個函數

import numpy as np
import matplotlib.pyplot as plt

def f1(x):
    y = np.exp(-x) * np.sin(2 * x)
    return y

def f2(x):
    y = np.exp(-x) * np.cos(3 * x)
    return y

x = np.linspace(0, 2 * np.pi, 100)
y1 = f1(x)
y2 = f2(x)

plt.plot(x, y1, '-', x, y2, '--')
plt.xlabel('x')
plt.ylabel('f(x) & g(x)')
plt.title('Plot of Two Functions')
plt.text(4, -0.2, 'Copyright@lin')

plt.show()

Laplace Transform

透過sympy進行Laplace transform

from sympy import *
from sympy.abc import t, s
laplace_transform(t, t, s)
inverse_laplace_transform(s**(-2), s, t)

註解: Heaviside(t)為單位步階函數\(u(t)\)

sympy質數

範例

[prime(i) for i in  range(1,10)]
isprime(10)
isprime(7)
prevprime(10)
nextprime(10)

輸出

>>> [prime(i) for i in  range(1,10)]
[2, 3, 5, 7, 11, 13, 17, 19, 23]
>>> isprime(10)
False
>>> isprime(7)
True
>>> prevprime(10)
7
>>> nextprime(10)
11

plt 字體異常

參見 Matplotlib 中文字體亂碼問題 · dw’s 小站

先檢查可用的字體列表

import matplotlib.font_manager
print([f.name for f in matplotlib.font_manager.fontManager.ttflist])

參見

Python