以下是一般任何數據分析, 我們推薦的標準套件。
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
sympy
和「sympy
化」¶我們先用一個簡單的例子來看符號型運算和傳統的數值型運算有什麼不同, 順便瞭解一下什麼叫 sympy
化。
今天我們主角是符號型運算的 sympy
(symbolic python)。符號型運算其實就是「數學課」用的運算方式, 例如 1/2 + 1/3 要告訴我們精確的 5/6, 而不是近似的數值解。
import sympy as sp
1/2 + 1/3
0.8333333333333333
在 sympy
裡, 最標準的方法是告訴電腦這是兩個有理數, 因此加出來也會是有理數。
sp.Rational(1,2) + sp.Rational(1,3)
但這看起來太煩, 我們要告訴電腦用 sympy
去算, 要用一個叫 sympify
(也就是 sympy
化) 的指令。
sp.sympify(1)/2 + sp.sympify(1)/3
但 sympify
也未免太長了, 所以有一個簡化到不能再簡化的版本 S
。
from sympy import S
S(1)/2 + S(1)/3
本來就是交給 sympy
算的我們自然是不用擔心。
sp.sqrt(2) * sp.sqrt(3)
考一下 sympy
三角函數程度如何。
π = sp.pi
π
sp.sin(π/3)
如果我們要看到近似的數值, 可以用 n
這個 method。
π.n(100)
$\pi$ 叫做超越數, 而很神奇的是世界上任何的數字, 包括你的手機號碼, 都會出現在 $\pi$ 中。文字也可以變成數字, 所以 $\pi$ 也會出你的名字等等!
你可以先試試你的生日, 比如 1215, 是不是有在 $\pi$ 裡面。
'1215' in str(π.n(100000))
True
電腦中的變數都要先給一個數值, 比如說
egg = 87
不然呼叫時會出現錯誤。但是數學中我們常常要用的變數, 像是
$$f(x) = x^2$$沒有指定值是多少, 甚至這樣就要給他微分或積分, 這該如何是好? 所以在 sympy
中就要有意識的告訴 Python, 這是我們要留給數學變數用的。定義方式有三種。
Symbol
¶假設我們要把 $x$ 定義為數學上的變數。
x = sp.Symbol('x')
x**2
是不是! Python 真的不會抱怨了! 但是為什麼這裡 $x$ 要寫兩次呢? 原因是前面的 x
是我們在電腦中要用的變數名稱, 後面的是顯示的名稱。在正常情況, 這兩者名稱當然是一樣的, 但有時我們也會設不一樣的名稱。
比如說下面這個例子。
x = sp.Symbol('🐶')
f = x**2
f
當然這麼耍寶的情況應該不多, 但有些情況很有用。比如我們想定 a 是顯示 $\alpha$...
a = sp.Symbol('alpha')
sp.sqrt(a)
symbols
一次令多個變數¶事實上我們很少用 Symbol
這個指令, 而是用更方便的版本, 一次令多個變數的 symbols
(沒錯, 變小寫加 s
!)
x, y, z = sp.symbols('x,y,z')
試用一下。
x**2 * y + 2*z
一定令多個希臘字母也可以!
a, b, c = sp.symbols('alpha,beta,gamma')
a*b + c**2 + 3*a*b*c
abc
子套件¶雖然偶而我們會想把變數名稱、顯示名稱選用不同的名字 (例如 a
這個變數, 顯示名稱是 $\alpha$。但很多情況我們會用相同的名稱, 那不管用 Symbol
或 symbols
都太麻煩了。
這裡我們介紹 abc
這個在 sympy
中的神秘子套件, 把
全部設成我們要當 sympy
用的變數。
from sympy.abc import x, y, z, t, a, b, c, alpha, beta, gamma
我們來試試看。
a*x + b*y + c*z + gamma*x**2
sympy
用的函數¶在 sympy
中要令一個函數 (事實上是定義一個式子) 很簡單。
f = x**2
f
subs
代值¶但要記得, 這不是一個 Python 的函式, 所以要計算 $f(87)$ 不能用
f(87)
要用 subs
代入。
f.subs(x, 87)
要是這函式有多個變數怎麼辦呢? 那就用一個串列。
g = x*y + 1
g.subs([(x, 94), (y, 87)])
我們甚至可以快速畫出這個函數的圖形。
sp.plot(f)
<sympy.plotting.plot.Plot at 0x7f0ccde30b10>
也可以指定 $x$ 的範圍。
sp.plot(f, (x, -5, 5))
<sympy.plotting.plot.Plot at 0x7f0cc5685710>
PlotGrid
同時顯示多個圖¶比較進階的是我們可以用 PlotGrid 同時顯式多個圖形。我們來用四大天王函數來試試。
from sympy.plotting import PlotGrid
p1 = sp.plot(x**2, show=False)
p2 = sp.plot(sp.exp(x), show=False)
p3 = sp.plot(sp.log(x), show=False)
p4 = sp.plot(sp.sin(x), show=False)
PlotGrid(2, 2, p1, p2, p3, p4)
<sympy.plotting.plot.PlotGrid at 0x7f0cc5617a90>
sympy
比求極限的功力¶sp.limit(x**2+1, x, 2)
哦, sympy
會耶! 但這麼簡單的實在沒有很令人感動。來個略略難一點的例子。
f = (x**2+2*x-3)/(x-1)
f
sp.limit(f, x, 1)
分段定義函數是可能的嗎? 比方說課本習題的一個例子:
$$f(x) = \begin{cases} 2x+1 & \mbox{當 } x\leq 2\\ 7-x & \mbox{當 } 2< x < 4\\ x & \mbox{當 } x\geq 4 \end{cases}$$f = sp.Piecewise((2*x+1, x<=2), (7-x, x<4), (x, True))
這個函數在 $x=4$ 是不連續的, 如果沒有不連續的點, sympy
要畫個分段定義函數還是強強的。但是有不連續的點, sympy
會把這點連起來! 為了看到正確的圖形, 我們用了分開畫的小技巧。
p1 = sp.plot(f, (x, 0, 3.99), show=False)
p2 = sp.plot(f, (x, 4, 6), show=False)
p1.extend(p2)
p1.show()
sp.limit(f, x, 5)
再來算關鍵點的部份!
sp.limit(f, x, 4, '-')
結果 sympy
算錯了!! 看右極限。
sp.limit(f, x, 4, '+')
這是正確的, 可惜單邊極限在 sympy
有 bug, 我們只好算更簡單的單邊極限。
sp.limit(1/x, x, 0, '+')
sp.limit(1/x, x, 0, '-')
順便說明一下, 在 sympy
要表示 $\infty$ 是用很可愛的 oo
。
f = (4*x**3-5)/(x**3-2*x+9487)
f
sp.limit(f, x, sp.oo)
到了看 sympy
微分功力的時刻了。
f = x**4 - 3*x**2 + 2*x - 5
f
sp.diff(f, x)
還可以啦, 畢竟很簡單。二階導數可以怎麼做呢?
sp.diff(f, x, x)
另外也可以這樣寫。
sp.diff(f, x, 2)
指數、對數函數也來試試。
sp.diff(sp.exp(x), x)
sp.diff(sp.log(x), x)
事實上三角函數也是可以的!
sp.diff(sp.sin(x), x)
sp.diff(sp.tan(x), x)
注意我們平常習慣看到的是 $\sec^2(x)$。
sp.diff(sp.sec(x), x)
看看好像很難的反三角函數的微分。
sp.diff(sp.asin(x), x)
sp.diff(sp.acos(x), x)
sp.diff(sp.atan(x), x)
這隱函數微分總會難倒 sympy
了吧? 我們來看個課本的例子。注意我們要給 sympy
$f(x,y)=0$ 這樣的 $f(x,y)$, 所以課本上的
我們要令成
$$y^3 + x^3y - 2xy^2 + 4x^2 - 4$$h = y**3 + x**3*y - 2*x*y**2 + 4*x**2 - 4
h
sp.idiff(h, y, x)
不得了, 這正確答案!! 我們當然也可以求課本上說這條式子定義出的曲線, 在 $x=1, y=0$ 的切線斜率。
sp.idiff(h, y, x).subs([(x,1), (y,0)])
完全正確! 好奇的話我們也可以畫一下這個式子定義出的曲線圖形。sympy
如果是
這種形式的隱函數, 要用這樣的方式輸入這個等式:
sp.Eq(f, a)
然後用 plot_implicit
畫出來。
sp.plot_implicit(sp.Eq(h, 0))
<sympy.plotting.plot.Plot at 0x7f0cc5492790>
solveset
找微分等於 0 的點¶我們知道一個函數的極值是發生在臨界點, 也就是微分等於 0 或是微分不存在的那些點。我們試試可愛的 sympy
可不可以幫我們算。
f = x**4 - 8*x**2 + 1
f
f1 = sp.diff(f, x)
sp.solveset(f1)
f2 = sp.diff(f, x, x)
f2.subs(x, -2)
f2.subs(x,0)
f2.subs(x,2)
由二階導數判別相對極值的方法, 可知函數在 $x=-2$ 時有相對極小值, 在 $x=0$ 時有相對極大值, 在 $x=2$ 時有相對極小值。
我們也學過, 函數的遞增遞減性只有在臨界點會產生變化。在我們這個例子, 我們可以判斷這幾個區間遞增遞減的情況:
$$(-\infty, -2), (-2,0), (0,2), (2, \infty)$$還記得就是區間中隨便找個代表點, 代入一階導數, 看正負號就可知函數在這個區間是遞增或是遞減的。
這個例子我們會知道函數是
遞減 -> 遞增 -> 遞減 -> 遞增
用 sympy
可以確認一下。
sp.is_decreasing(f, sp.Interval(-sp.oo, -2))
True
sp.is_increasing(f, sp.Interval(-2, 0))
True
sp.is_decreasing(f, sp.Interval(0, 2))
True
sp.is_increasing(f, sp.Interval(2, sp.oo))
True
使用時要小心一個狀況, 假設我們的函數是這個:
$$f(x) = x^3 + x$$f = x**3 + x
f
f1 = sp.diff(f, x)
f1
sp.solveset(f1, x)
等一下, 這在實數上是沒有解的, 所以無局部極大或極小。雖然我們這樣也看得出來, 但能不能讓 sympy
告訴我們實數上是無解的呢?
sp.solveset(f1, x, domain=S.Reals)
singularities
找出微分不存在的點¶微分不存在的點呢? 我們可以用 singularities
來確認是不是有微分不存在的點。首先先看我們標準例子
from sympy.calculus import singularities
f = sp.Abs(x)
f
f1 = sp.diff(f, x)
singularities(f1, x)
f = 1/x + x
f1 = sp.diff(f, x)
singularities(f1, x)
終於到了我們積分的時刻。首先, 我們先來個假的積分, 我是說, 不定積分。
f = sp.sqrt(31*x + 2)
f
sp.integrate(f, x)
答案是正確的, 但是沒有加 C 要扣分...
再來看定積分。
f = x/(sp.sqrt(x**2+4))
f
sp.integrate(f, (x,0,2))
事實上這是 sympy
, 你要看一般從 a 到 b 的積分值也可以。
sp.integrate(f, (x,a,b))
習題抽幾題考考 sympy
。
sp.integrate(sp.sqrt(sp.log(x))/x, x)
sp.integrate(x*sp.exp(-x**2), (x, 0, 1))
還記得微積分基本定理有一個形式, 告訴我們微分和積分基本上是反運算。寫成數學式子是這樣:
$$\dfrac{d}{dx} \int_a^x f(t) \, dt = f(x)$$這寫起來就好可怕的東西, sympy
真的會嗎? 我們來試試。
f = t**3/(t**2+1)
f
sp.diff(sp.integrate(f, (t, a, x)),x)
檢查一下, 這是正確答案。自然我們可以叫 sympy
幫我們化簡整理。
sp.simplify(sp.diff(sp.integrate(f, (t, a, x)),x))
至此都是「簡單的」積分。那需要用到積分大魔王「部分積分」技巧的, sympy
還會做嗎?
sp.integrate(x*sp.exp(x), x)
sp.integrate(sp.log(x), x)
好像都難不倒 sympy
, 那最後我們來看看暇積分吧。
sp.integrate(1/x**2, (x, 1, sp.oo))
sp.integrate(2*x*sp.exp(-x**2), (x, -sp.oo, sp.oo))
最後我們來看超級酷炫, 把可以一直微分的函數轉成多項式的技術: 泰勒級數!
一個函數 $f(x)$ 在 $x=a$ 的泰勒展開式是長這個樣子:
$$f(x) = f(a) + f'(a) (x-a) + \dfrac{f''(a)}{2!} (x-a)^2 + \cdots + \dfrac{f^{(n)}(a)}{n!} (x-a)^n + \cdots$$這裡順便說一下, 我們差不多從微積分的第一天開始, 就灌輸大家一個概念: 一個函數 $f(x)$ 如果在 $x=a$ 那點的切線存在的話 (也就是可微), 那這個函數的圖形 $y=f(x)$ 在 $x=a$ 附近和那條切線長得很像。也就是說當 $x$ 和 $a$ 很接近的時候, 我們有:
$$f(x) \simeq f(a) + f'(a) (x-a)$$耶, 這就是泰勒展開式的第一項啊! 如果我們再用二階導數再調整一點就會更像, 三階再調一點, ... 這就得到泰勒級數!
說得這麼多, 結果 sympy
要求泰勒展開式很容易。我們來看看我們四大天王中的
在 $x=0$ 的泰勒展開式。
f = sp.exp(x)
sp.series(f, n=10)
最後一項是說真正的 $e^x$ 和這個泰勒展開式差距最多是 $x^{10}$ 的某個固定的常數倍。
removeO
去掉誤差項¶計算的時候我們希望不要有這個誤差估計項, 那要怎麼做呢?
f_taylor = sp.series(f, n=10).removeO()
f_taylor
把一個函數變成可愛的多項式, 就真的可以算了! 比如說, 這裡我們可以開心的算一下 $e = e^1$ 的值。
e = f_taylor.subs(x,1)
e
等等, 這個數字好像沒什麼意義。不過我們前面說過換成數值解的方法。
e.n(100)
好奇的話我們可以比較一下正確答案。
sp.exp(1).n(100)
我們來用互動感受一下泰展多項式越來越像原來函數的感覺。
from ipywidgets import interact
def taylor(n=10):
f = sp.exp(x)
egg = sp.series(f, n=n)
f_taylor = egg.removeO()
p1 = sp.plot(f, (x, -5, 5), line_color='blue', show=False)
p2 = sp.plot(f_taylor, (x, -5, 5), line_color='red', show=False)
p1.extend(p2)
p1.show()
taylor(n=10)
interact(taylor, n=(1,50))
interactive(children=(IntSlider(value=10, description='n', max=50, min=1), Output()), _dom_classes=('widget-in…
<function __main__.taylor>
我們知道
$$\tan(\pi/4) = 1$$sp.tan(π/4)
也就是說
$$\arctan(1) = \pi/4$$於是我們求出 $\arctan(x)$ 的泰勒多項式, 計算 $\arctan(1)$ 的值, 乘以 4 就是 $\pi$!。
f = sp.atan(x)
n = 100
f_taylor = sp.series(f, n=n).removeO()
pi = 4*f_taylor.subs(x, 1)
pi
pi.n(100)
等等, 這小數點第二位就錯了, 這實在太遜了一點! 沒有錯, 其實這是收斂超慢計算 $\pi$ 的方法。比較好的方法像是數學神人拉馬努金 (Ramanujan) 我們根本猜不到他怎麼想出來, 和 $\pi$ 有關的等式:
$$\dfrac{4}{\pi} = \dfrac{1}{882} \sum_{n=0}^\infty \dfrac{(-1)^n(4n)!}{(4^nn!)^4} \dfrac{1123 + 21460n}{882^{2n}}$$有興趣可以用 sympy
算算看, 會發現比我們上面用的方法好太多!