1.極坐標系統（Polar coordinate system）
#（1）注意：sympy對於極坐標的支援有限：
#所以若是要畫出極坐標圖，就只能用直角坐標的x,y來畫圖
#極坐標與直角坐標的轉換公式：
# (r*cos(theta), r*sin(theta)) for x and y coordinates

#(2).極坐標系統： p(r,phi)
#極坐標與直角坐標的轉換：x=r*cos(theda)，y=r*sin(theda)

範例4-1：極坐標系統（Polar coordinate system）

2.範例4-1：極坐標系統（Polar coordinate system）
#sympy沒有提供polar的網格與運算，只能用plot_parametric畫(x,y)
#cmath + matplotlib模組，有提供ploar網格，可以直接畫極坐標圖

from sympy import *
#(1).若用symbols+plot()：則畫出正弦圖種類圖,不是圓
theda = symbols('theda')
f = '5*cos(theda) + 5*sin(theda)'
plot(f,(theda,0,pi))
#print('sin(pi/6)=', f.subs(theda,pi/6))

#(2).若用plot_parametric：則可畫出極坐標圖，是圓
theda = symbols('theda')
#plot_parametric(x軸值, y軸值,(變數範圍,-5,5))
#用複數的概念，用plot_parametric(x, yi)畫圖
plot_parametric(cos(theda), sin(theda), (theda, -5, 5))

成果圖示：
成果圖示：

程式碼內容

範例4-2：複數系統cmath的應用：三種坐標轉換（複數坐標，極坐標，直角坐標）

3.範例4-2：複數系統cmath的應用：三種坐標轉換（複數坐標，極坐標，直角坐標）
#(1).Python的兩個內建數學模組: math, cmath
# math 模組：適合處理非向量運算，且運算效能較好，但math 模組只能處理實數運算
#cmath模組： 專門用來處理複數運算
#(2).cmath 模組的函數 :
#cmath模組的運算元是複數, 常數增加了複數的 infj 與 nanj, 方法則大部分與 math 相同, 不過沒有階乘函數 (因為複數沒有階乘),
#新增三個複數座標轉換的方法（在卡式座標與極座標之間轉換）:
#polar(x) 傳回複數 x 的極座標表示法元組 (r, p), r=長度, p=角度
#rect(r, p) 傳回極座標 (r, p) 的複數
#phase(x) 傳回複數 x 的弧度 (radian)
import cmath, math

print('全部cmath的指令 = ', dir(cmath))

#已知複數，傳回弧度phase（theda）
d1 = cmath.phase(5+5j)
print('已經複數，傳回弧度theda=', d1)
print('再轉成角度=', math.degrees(d1) )

#已知複數，轉成極座標:polar(a+bi) =》 p(r,theda)
print('複數轉極座標(r,theda)=', cmath.polar(5+5j))

#已知極座標，轉成直角坐標：p(r,theda) => rect(x,y)
print('極座標轉直角坐標：p(r,theda) => rect(x,y)=',cmath.rect(7.0710678118654755, 0.7853981633974483))

#求直角坐標(1,1)的複數
print('求直角坐標(1,1)的複數=',complex(1.0,1.0))
#求直角坐標(1,1)的弧度
print('求直角坐標(1,1)的弧度=', cmath.phase(complex(1.0, 1.0)))
#求直角坐標(1,1)的角度
print('求直角坐標(1,1)的角度=', math.degrees(cmath.phase(complex(1.0, 1.0))))
#求直角坐標(1,1)的極坐標
print('求直角坐標(1,1)的極坐標p(r,theda)=',cmath.polar(complex(1.0, 1.0)))

#cmath.polar() 與 cmath.rect() 在坐標轉換時，非常實用

成果圖示：

程式碼內容

範例4-3：虛數的物理意義，虛數相乘的物理意義

4.範例4-3：虛數的物理意義，虛數相乘的物理意義
#複數：a+bj
#注意，它使用j，不是i

import cmath
cpx = 2 + 3j
print('2+3j=', cpx)

#虛數定義：i² = -1
print('j*j=',1j*1j)
print('cmath.sqrt(-1) =',cmath.sqrt(-1))

#i 的本質是單位週期結構最基本形式，它並不是一個數，確切地說就是一個旋轉量
#說明虛數的物理意義“
#(1).實數1向左旋轉90度後，會變成i
#(2). i再旋轉90度，變成-1
#(3).這個旋轉90度，就是i
#(4).i 就是意味著逆時針旋轉90°，-i 就是順時針旋轉90°
#(5).若向左旋轉45度，就變成1+1i
#(6).虛數相乘的物理意義：改變旋轉角度
#例如：箭頭是 3 + 4i，若向左旋轉45°（1+1j)
#虛數相乘 = （3+4j)*(1+1j)) = -1+7j
#(7)虛數相乘的極坐標表示 = (a+bi)(c+di)＝r1*r2*(cos(α+β) + isin(α+β))
#(8)虛數相加的物理意義：合成向量

#實部：real，和虛部：imag
print('實數=', cpx.real)
print('實數=', cpx.imag )

#共軛數：conjugate
#複數是複平面上的一個向量，而共軛複數關於實軸對稱（上下對稱）
print('共軛數=', cpx.conjugate())

#求複數的絕對值(a^2+b^2)
print('複數的絕對值=', abs(cpx))

#把實數取複數
print('把實數2取複數=', complex(2))
print('把實數（2,3)取複數=', complex(2,3))
print('設定real=2,imag=3，取複數=', complex(real=2,imag=3))

成果圖示：

程式碼內容

範例4-4：在極坐標上畫出一個點，畫出一條線

5.範例4-4：在極坐標上畫出一個點，畫出一條線
import numpy as np
import matplotlib.pyplot as plt

#畫出極坐標圓形圖，要設定axes的projection=‘polar'
ax = plt.axes(projection='polar')

# 設定標題
plt.title('Circle of polar:')

#畫出一個點
#公式ax.plot(r,theda,'mark')
ax.plot(np.pi, 0.5,'o')

#根據三個點，畫出線條
theta = [0, 1, 2.5]
r = [0, 0.5, 1]
ax.plot(theta, r)
#若要加上點，就要加上'o'
#ax.plot(theta, r,'o')

#存檔
#plt.savefig("out.png")

成果圖示：

程式碼內容

範例4-5：畫出極坐標的曲線方程式

6.範例4-5：畫出極坐標的曲線方程式
import cmath
import matplotlib.pyplot as plt
import numpy as np

ax = plt.axes(projection='polar')
#設定曲線方程式p(r,theda)，其中r是變數，逐漸變大
r = np.arange(0, 12, 0.1)
theda = r * cmath.pi
#ax.plot(r,theda)
ax.plot(r,theda,'o',c='red', alpha=0.35)

成果圖示：

程式碼內容

範例4-6：用sympy + plot3d_parametric_line畫3D曲線

7.範例4-6：用sympy + plot3d_parametric_line畫3D曲線
from sympy import *
u = symbols('u')
#原理：本來若只畫plot_parametric_line(x,y)，就只是一個極坐標圓
#若是畫出三個參數，可用plot3d_parametric_line(x,y,z),其中z=u，就會把極坐標圓拉開，拉長
#語法：plot3d_parametric_line(x, y, z, (變數, -5, 5))
plotting.plot3d_parametric_line(cos(u), sin(u), u, (u,-5,5))

成果圖示：

程式碼內容

範例4-7：用sympy + plot3d_parametric_line畫3D曲面(surface)

8.範例4-7：用sympy + plot3d_parametric_line畫3D曲面(surface)
from sympy import *
u,v = symbols('u v')

#plotting.plot3d_parametric_surface(cos(u+v), sin(u-v),u,(u,-5,5),(v,-5,5))
plotting.plot3d_parametric_surface(v**(2), u**(2), u, (u,-5,5),(v,-5,5))

成果圖示：

程式碼內容

1.一階微分，一階導數

1.一階微分，一階導數
牛頓符號 = f'(x)
萊布尼茨 符號 = df(x)/dx = dy/dx
df(x)/dx = dy/dx = lim,x→0 (f(x)-f(a))/(x-a)

範例：求 y=1/x = x^-1的微分
dy/dx = -2*x^-2

範例：求 y= 3*x^2 + 2x + 5的微分
dy/dx = 6*x1 +2

一階微分的物理意義 = 代表曲線的斜率
一階微分的應用：微分大於0，代表斜率為正。微分小於0，代表斜率為負
一階微分的應用：斜率=0處（由正轉負，或由負轉正），就是曲線局部有極值（極大值，或極小值）

範例5-1：求微分，常見的幾類微積分函數基本公式

2.範例5-1：求微分，常見的幾類微積分函數基本公式
from sympy import *
x=symbols('x')
f = Function('f')(x)
#一階導數(對x) = 一階微分(對x) = diff(函數, x)
f = 1/x
print('求 y=1/x的微分=', diff(f, x))

f = 3*x**2 + 2*x + 5
print('求 y=3*x^2 + 2*x + 5的微分=', diff(f, x))

#常見的幾類微積分基本公式
f = sin(x)
print('求sin(x)的微分=', diff(f, x))
f = cos(x)
print('求cos(x)的微分=', diff(f, x))

#標準指數e的微分
f = E**x
print('求E**x的微分=', diff(f, x))
f = exp(x)
print('求exp(x)的微分=', diff(f, x))

#自然對數e的微分
#自然對數 = ln()
#10為底的對數 = log()
f = ln(x)
print('求自然對數ln(x)的微分=', diff(f, x))

#多項式的微分
#如果方程式有未知參數n,就必須要加上''
f = 'x**n'
print('求x**n的微分=', diff(f, x))
f = x**5
print('求x**5的微分=', diff(f, x))

成果圖示：
程式碼內容

範例5-2：求在某點位置的微分值（斜率）

3.範例5-2：求在某點位置的微分值（斜率）
from sympy import *

x = symbols('x')
f = Function('f')(x)

f = 1/x
#帶入指定行數值的方法有兩種：.subs({x:2})。或是.subs(x,2)
print('1/x在x=2的微分值=', diff(f, x).subs({x:2}))
print('1/x在x=2的微分值=', diff(f, x).subs(x,2))

f = 3*x**2 + 2*x + 5
print('求 y=3*x^2+2*x+5，在x=2的微分=', diff(f, x).subs(x,2))

#常見的幾類微積分基本公式
f = sin(x)
print('求sin(x)在pi/6的微分值=', diff(f, x).subs(x,pi/6))
f = cos(x)
print('求cos(x)在pi/6的微分值=', diff(f, x).subs(x,pi/6))

#標準指數e的微分
f = E**x
print('求E**x在x=2的微分值=', diff(f, x).subs(x,2))
f = exp(x)
print('求exp(x)在x=2的微分值=', diff(f, x).subs(x,2))

#自然對數e的微分
#自然對數 = ln()
#10為底的對數 = log()
f = ln(x)
print('求自然對數ln(x)在x=2的微分值=', diff(f, x).subs(x,2))

成果圖示：
程式碼內容

4.二階微分，二階導數

4.二階微分，二階導數
牛頓符號 = f''(x)
萊布尼茨 符號 = (d/dx) df(x)/dx = d^2y/dx^2

範例：求 y= 3*x^2 + 2x + 5的二次微分
dy/dx = 6*x1 +2
f''(x) = 6

二階微分的物理意義 = 斜率變化率（可以用找出反曲點）
二階微分的應用：二階微分為正數,表示是向上開口形狀的曲線。
二階微分為負數，表示是向下開口形狀的曲線

範例5-3：求二階導數，判別反曲點，判別線型凹凸

5.範例5-3：求二階導數，判別反曲點，判別線型凹凸
#求 y= 3*x^2 + 2x + 5的二次微分
from sympy import *
x = symbols('x')
f = Function('y')(x)

f = 3*x**2 + 2*x + 5

#高階的微分可以使用 diff(func, var, n)
print('求 y= 3*x^2 + 2x + 5的二次微分=', diff(f, x, 2))

f = -x/(x**2+4)
#求出這個曲線的反曲點，與凹凸
df2f2 = diff(f, x, 2)
print('求 y= 3*x^2 + 2x + 5的二次微分=', diff(f, x, 2))
#二階導數 = 0處，就是反曲點
print('求二階導數 = 0的解，',solve(df2f2))
#若 x <0-2*sqrt(3))，由二階導數的值正負，判別曲面上下
if diff(f, x, 2).subs(x,-2*sqrt(3)-1) <0:
print('x <-2*sqrt(3)，二階導數的值為負，曲面向下')
else:
print('x <-2*sqrt(3)，二階導數的值為正，曲面向上')

#若 -2*sqrt(3) if diff(f, x, 2).subs(x,-1) <0:
print('-2*sqrt(3) else:
print('-2*sqrt(3)
#若 0 if diff(f, x, 2).subs(x,1) <0:
print(' 0 else:
print(' 0
#若 2*sqrt(3) < x，<由二階導數的值正負，判別曲面上下
if diff(f, x, 2).subs(x,5) <0:
print(' 2*sqrt(3) < x，二階導數的值為負，曲面向下')
else:
print(' 2*sqrt(3) < x，二階導數的值為正，曲面向上')

#畫出f2的線性
plot(f)
成果圖示：
成果圖示：
程式碼內容

範例5-4：求一階導數，判別局部極大值，極小值的位置

6.範例5-4：求一階導數，判別局部極大值，極小值的位置
from sympy import *
x = symbols('x')
y = Function('y')(x)
y = sqrt(3)*cos(x) + sin(x) + 1
#求f的局部極值（極大值，極小值）
dif1y = diff(y, x, 1)
print('一階導數函數 = ', dif1y)
#求解一階導數=0的位置
print('求解一階導數=0的解 = ', solve(dif1y,x))
#兩個解，x0,x1
x0 = solve(dif1y,x)[0]
x1 = solve(dif1y,x)[1]
print(x0,x1)

#判別斜率正負
if dif1y.subs(x,x0-0.1) < 0 and dif1y.subs(x,x0+0.1) > 0:
print('x=',x0,'的位置是局部最小值')
elif dif1y.subs(x,x0-0.1) >0 and dif1y.subs(x,x0+0.1) <0:
print('x=',x0,'的位置是局部最大值')

if dif1y.subs(x,x1-0.1) < 0 and dif1y.subs(x,x1+0.1) > 0:
print('x=',x1,'的位置是局部最小值')
elif dif1y.subs(x,x1-0.1) >0 and dif1y.subs(x,x1+0.1) <0:
print('x=',x1,'的位置是局部最大值')

#畫出y的線性
plot(y, line_color='red')

成果圖示：
成果圖示：
程式碼內容

1.積分定義

1.積分定義：
函數 y = f(x)
不定積分 = F(x) = ∫f(x)dx = 算式 = 方程式（沒有給定上限b, 下限a)
定積分 = F(x) = ∫f(x)dx （給定上限b, 下限a)

範例：y = 2x^2 + 3x
求y的積分（0,2) = F(2) - F(0)
(1).不定積分 F(x) = 2/3*x^3 +3/2*x^2 + C
這個常數C有無限多種可能，但是求定積分時，F(2)-F(0) 會把C刪除
(2).定積分時，F(2)-F(0) = ...

(3).積分的物理意義 =
第一種物理意義：定積分乃是反推原函數的單位時間(x）的變化率（F1→ F2）
（原函數 = 積分函數的變化率）
（因為：積分函數的導數就是原函數，原函數的反導函數就是積分）
（積分的導函數 = 變化率 = 例如：速度s/t，或是金融的邊際利潤率）
（例如F是位移，f=速度=F/t =單位時間內的距離改變 = 速度）

第二種物理意義：定積分是(0~2)的函數下方的面積

第三種物理意義：積分代表機率
（若f(x)dx的積分綜合　= 1，那麼定積分（x=a~b）間的面積總和積分，就代表機率
（f(x)代表機率密度函數）

(4).積分（反導函數）有無限多組的可能解
F(x) = 2/3*x^3 +3/2*x^2 + C1
F(x) = 2/3*x^3 +3/2*x^2 + C2
F(x) = 2/3*x^3 +3/2*x^2 + C3
...

2.比較：積分，函數，微分，二階微分的關係：（導數 vs 反導數）

2.比較：積分，函數，微分，二階微分的關係：（導數 vs 反導數）
(1).積分的導數 = 函數
函數的導數 = 一階微分
一階微分的導數 = 二階微分

(2).積分的導數 = 函數
二階微分的反導數 = 一階微分
一階微分的反導數 = 函數
函數的反導數 = 積分

3.sympy的積分語法：integrate

3.sympy的積分語法：integrate
（1）不定積分：integrate(f, x)
returns the indefinite integral ∫fdx∫fdx

（2）定積分：integrate(f, (x, a, b))
returns the definite integral ∫bafdx∫abfdx

範例6-1：求函數的不定積分，定積分

4.範例6-1：求函數的不定積分，定積分
from sympy import *
x = symbols('x')
f = Function('f')(x)
#計算y = 2x^2 + 3x的不定積分
#語法：integrate(f1,x)
f = 2*x**2 + 3*x
print('求y = 2x^2 + 3x的不定積分=', integrate(f,x))
#注意：最後的積分還要自己手動去加上任意常數C

#計算y = 2x^2 + 3x的定積分（0～2）
#語法：integrate(f1,(x,0,2))
f = 2*x**2 + 3*x
print('求y = 2x^2 + 3x的不定積分=', integrate(f,(x,0,2)))

#輸出標準數學符號
from sympy.interactive import printing
printing.init_printing(use_latex=True)
display(f)

成果圖示：
程式碼內容

範例6-2：比較sympy與Maxima軟體在計算積分的準確率

5.範例6-2：比較sympy與Maxima軟體在計算積分的準確率
from sympy import *
from sympy.interactive import printing
printing.init_printing(use_latex=True)

x = symbols('x')
y1 = Function('y')(x)
#使用軟體或函數庫來計算微分都正確，但是若是計算積分，很多函數會算不出來
#微積分常用軟體Maxima在做積分時，經常會出現看不懂的函數結果，或是算不出答案

#(1).簡單積分
y1 = x/sqrt(1+x**2)
print('1.簡單積分=',integrate(y1, x))
display(y1)
#注意：最後的積分還要自己手動去加上任意常數C

#(2).簡單積分
y1 = 1/x
print('2.簡單積分=', integrate(y1, x))
display(y1)
#注意：最後的積分還要自己手動去加上任意常數C

#(3).sympy與Maxima軟體的結果一樣複雜
y1 = sqrt(1+x**2)
print('3.sympy與Maxima軟體的結果一樣複雜=', integrate(y1, x))
display(y1)
#結果：與Maxima軟體一樣，都會出現asinh(x)這個少見的函數

#(4).sympy可以積分出來，但Maxima軟體無法正確積分
y1 = sqrt(1+x**4)
print('4.sympy可以積分出來，但Maxima軟體無法正確積分=\n', integrate(y1, x))
display(y1)
#結果：若用Maxima軟體計算結果無法顯示答案，但若用sympy則可以顯示出答案

#計算0～1之間的面積（積分）
print('5.計算0～1之間的面積=\n', integrate(y1,(x,0,1)))
display(y1)

成果圖示：
程式碼內容

範例6-3：三角函數的積分

6.範例6-3：三角函數的積分
from sympy import *
u = symbols('u')
f = Function('f')(u)

#1.sin(u)的積分 = -cos(u)
f = sin(u)
print('1.sin(u)的積分 = -cos(u)=', integrate(f,u))
#注意：最後的積分還要自己手動去加上任意常數C

#2.cos(u)的積分 = -sin(u)
f = cos(u)
print('2.cos(u)的積分 = sin(u)=', integrate(f,u))
#注意：最後的積分還要自己手動去加上任意常數C

#3.tan(u)的積分
f = tan(u)
print('3.tan(u)的積分=', integrate(f,u))
#注意：最後的積分還要自己手動去加上任意常數C

#4.sec(u)的積分
f = sec(u)
print('4.sec(u)的積分=', integrate(f,u))
#注意：最後的積分還要自己手動去加上任意常數C

成果圖示：
程式碼內容

1.常微分方程式(ordinary differential equation，簡稱ODE）

1.常微分方程式(ordinary differential equation，簡稱ODE）C
（1）微分方程式是未知函數只含有一個自變數的微分方程式

（2）方程式：
方程式 = 函數變數的等式，就可以求解
一元一次方程式：3y = 5
一元二次方程式：x**2 + 2x = 6
二元一次方程式：2x - y = 5

（3）代數方程式：
代數方程式：方程式的係數，解，都是『數』

（4）微分方程式：
微分方程式的解，不是個『數』，而是個『函數』

一階微分方程式：含有一階微分的方程式
例如：∂y/∂x = x

二階微分方程式：含有二階微分的方程式
例如：∂/∂x(∂y/∂x) = x
C
（5）微分方程式的解是個函數，而且有無限多組的可能解
例如：∂y/∂x = 1
這個解，是y=x，是個函數

例如：∂y/∂x = x
這個解，是y=(x**2)/2，是個函數，這個解不是唯一。
因為以下都可以是解答。
y=(x**2)/2 + 1
y=(x**2)/2 + 2
y=(x**2)/2 + 3
......

這個與積分（反導函數）有無限多組的可能解，是一樣的

（6）求微分方程式的解，就與算積分（反導函數）是一樣的
積分方程式：y = F(x) = ∫f(x)dx
微分方程式：∂y/∂x = f(x)
這兩個是一樣的，只是要把微分方程式→ 推算反導函數，就變成積分方程式了

所以，它們兩個都是有無限多組的可能解

2.求解一階微分方程式的唯一解，必須給定一個邊界條件

2.求解一階微分方程式的唯一解，必須給定一個邊界條件
例如：
距離 = y = F(x) = ∫f(x)dx
速度 = ∂y/∂x = f(x)

例如：目前已經微分（速度）=80m/s，請問車子的位置（距離）
結果：這個解答，有無限多組的可能，
除非：給定邊界條件，才能知道精確解答位置（唯一解)
邊界條件：y(a) = b

例如：微分方程式：y'(x) = ∂y/∂x = x
邊界條件是： y(1) = 1
求解？
先把微分方程式→積分→ 結果： y = ∫f(x)dx = (x**2)/2 + C
代入邊界條件：y(1) = 1/2+C = 1
知道，C = 1/2
求出微分方程式的唯一解答(函數）→ y = (x**2)/2 + 1/2

3.求解二階微分方程式的唯一解，必須給定兩個邊界條件

3.求解二階微分方程式的唯一解，必須給定兩個邊界條件
例如：加速度，速度，距離的關係式：
加速度 y''(t) = g = ∂/∂t(∂y/∂t) = 9.8
第一個邊界條件：起始速度 = y'(0) = -2
第二個邊界條件：起始位置 = y(0) = 10

求解步驟：
A.先把二階導數積分：y'(t) = ∫-9.8dt = -9.8x + C1（是個函數）
B.給定第一個邊界條件：y'(0) = -2
得到 C1 = -2
C.再把一階導數函數積分 y = ∫(-9.8t -2)dt = -4.9t**2 - 2t + C2
D.給定第二個邊界條件： y(0) = 10
得到 C2 = 10
E.得到微分方程式的唯一解 → y = -4.9t**2 -2x + 10

4.sympy的微分函數寫法有兩種

4.sympy的微分函數寫法有兩種：
(1).方法一：y = Function('y')(x)
from sympy import *
x = symbols('x')
y = Function('y')(x)
ans = dsolve(ode, y, ics={y.subs(x,1):1})
print('ODE求解 ： ', ans.args[0],' = ', ans.args[1])

(2).方法二：y = Function('y')
from sympy import *
x = symbols('x')
y = Function('y')
ode = diff(y(x),x) -x**3 -1
ans = dsolve(ode,y(x),ics={y(0):3})
print('求出解：',ans.args[0],' = ', ans.args[1])

(3).比較：這兩個方法最大的不同，是在邊界條件的設定
ics={y.subs(x,1):1}
或是
ics={y(0):3

範例7-1：微分方程式ODE求解，y'= ∂y/∂x = x

5.範例7-1：微分方程式ODE求解，y'= ∂y/∂x = x
#微分方程式：y'(x) = ∂y/∂x = x
#邊界條件是： y(1) = 1
#原始方程式： y(x) = (x**2)/2 + 1/2
#降階，微分方程式ODE, y'(x) = x
#語法：dsolve(eq, f(x), ics={y(1):1})

from sympy import *
#1.第一種寫法
x = symbols('x')
y = Function('y')(x)
#注意：ode是個方程式，不是函數
ode = Eq(y.diff(x), x)

print('1.第一種寫法：')
print('沒有給定邊界條件是，C1無限多組可能')
#from ode to find the 原始函數 y(x)
ans = dsolve(ode, y)
print('ODE解是個函數=',ans)
print('ODE求解 ： ', ans.args[0],' = ', ans.args[1])
#邊界條件作為字典傳遞給 dsolve，通過 ics 命名參數
ans = dsolve(ode, y, ics={y.subs(x,1):1})
print('2.給定邊界條件是： y(1) = 1')
print('ODE解是個函數=',ans)
print('ODE求解 ： ', ans.args[0],' = ', ans.args[1])

#2.第二種寫法：
x = symbols('x',real = True) #實數
y = symbols('y',cls=Function)(x) #函數 y=x**2.....
#注意：ode是個方程式，不是函數
ode = Eq(y.diff(x), x)
#from ode to find the 原始函數 y(x)
ans = dsolve(ode, y, ics={y.subs(x,1):1})
print('\n3.第二種寫法：給定邊界條件是： y(1) = 1')
print('ODE求解，', ans.args[0],'=', ans.args[1])

#顯示標準數學符號
from sympy.interactive import printing
printing.init_printing(use_latex=True)
display(y)

#畫出y(x)
f = ans.args[1]
plot(f)

成果圖示：
第一種微分寫法程式碼
第二種微分寫法程式碼

範例7-2：微分方程式ODE求解，y'(x) = x**3+1， 邊界條件：y(0)=3

5.範例7-2：微分方程式ODE求解，y'(x) = x**3+1， 邊界條件：y(0)=3
#微分方程式：y'(x) = x**3+1
#邊界條件是： y(0)=3
#原始方程式： y(x) = x**4 + x + 3
#降階，微分方程式ODE, y'(x) = x**3+1
#語法：dsolve(eq, f(x), ics={y(0):3})
from sympy import *
x = symbols('x')
y = Function('y')(x)
#注意：ode是個方程式，不是函數
ode = Eq(y.diff(x), x**3 +1)
#from ode to find the 原始函數 y(x)
ans = dsolve(ode,y,ics={y.subs(x,0):3})
print('求出解：',ans.args[0],' = ', ans.args[1])

#畫出y(x)
f = ans.args[1]
#plot(f)
plot(f,(x,-4,4))

成果圖示：
成果圖示：
成果圖示：
第一種微分寫法程式碼
第二種微分寫法程式碼

範例7-3：微分方程式ODE求解，y'(x)-2x/y=0， 邊界條件：y(1)=2

6.範例7-3：微分方程式ODE求解，y'(x)-2x/y=0， 邊界條件：y(1)=2
#微分方程式：y'(x)-2x/y=0
#邊界條件是： y(1)=2
#原始方程式： y(x) = sqrt(2*x**2 + 2)
#降階，微分方程式ODE, y'(x)-2x/y=0
#語法：dsolve(eq, f(x), ics={y(1):2})
from sympy import *
x = symbols('x')
y = Function('y')
#注意：方程式裡面的y,都要寫成y(x)
ode = diff(y(x),x) - (2*x)/y(x)
#from ode to find the 原始函數 y(x)
ans = dsolve(ode, y(x), ics={y(1):2})
print('求出解：',ans.args[0],' = ', ans.args[1])

#畫出y(x)
f = ans.args[1]
plot(f)
成果圖示：
成果圖示：
第一種微分寫法程式碼
第二種微分寫法程式碼

範例7-4：微分方程式ODE求解，y'(x) = 3x**2(y+2)， 邊界條件：y(0) = 8

7.範例7-4：微分方程式ODE求解，y'(x) = 3x**2(y+2)， 邊界條件：y(0) = 8
#微分方程式：y'(x) = 3x**2(y+2)
#邊界條件是： y(0) = 8
#原始方程式： y(x) = 10*exp(x**3) - 2
#降階，微分方程式ODE, y'(x) = 3x**2(y+2)
#語法：dsolve(eq, f(x), ics={y(0):8})
from sympy import *
x = symbols('x')
#我的目的是要找出原始函數：y(x)
y = Function('y')
ode = diff(y(x),x) -3*x**2*(y(x) + 2)
ans = dsolve(ode, y(x), ics={y(0):8})
print('求解ode的原始函數：', ans.args[0],' = ', ans.args[1])

#畫出圖（但是似乎smbpy無法畫出正確圖）
f = ans.args[1]
plot(f)

成果圖示：
第一種微分寫法程式碼
第二種微分寫法程式碼

範例7-5：試求在 xy 平面中通過點(1,1)且斜率為-y/x之曲線

8.範例7-5：試求在 xy 平面中通過點(1,1)且斜率為-y/x之曲線
#題目的條件1：斜率 y'(x) = -y/x
#題目的條件2：邊界條件，通過 y(1) =1
#微分方程式：y'(x) = -y/x
#邊界條件是： y(1) = 1
#原始方程式： y(x) = 10*exp(x**3) - 2
#降階，微分方程式ODE, y'(x) = -y/x
#語法：dsolve(eq, f(x), ics={y(1):1})
from sympy import *
x = symbols('x')
y = Function('y')(x)
#注意：ode是個方程式，不是函數
ode = Eq(y.diff(x), -y/x)
ans = dsolve(ode, y, ics={y.subs(x,1):1})
print('ode求出解：',ans.args[0],'=',ans.args[1])

#畫出圖
f = ans.args[1]
#sympy若要畫出紅色線條: line_color = 'red'
plot(f,line_color='red')

成果圖示：
成果圖示：
第一種微分寫法程式碼
第二種微分寫法程式碼

範例8-1：二階常係數齊性 ODE 的解法(一)--相異實根

1.範例8-1：二階常係數齊性 ODE 的解法(一)--相異實根
# ode ： y''(x) + y'(x) -6y = 0
#邊界條件或起始條件：y(0) =1, y'(0) =2
#語法：dsolve(eq, f(x), ics={y(0):1,y(x).diff(x).subs(x,0):2})
#所解出的原方程式： y(x) = exp(2*x)

from sympy import *
x = symbols('x')
y = Function('y')(x)
#注意：ode是個方程式，不是函數
ode =Eq(y.diff(x,2) + y.diff(x) - 6*y, 0)
#沒有考慮邊界條件下的通解：
ans = dsolve(ode, y)

#考慮邊界條件下的的特解(Particular Solution)：
#ics語法((1)：ics={y.subs(x,0):1,y.diff(x).subs(x,0):2})
#ics語法(2)：ics={f(x0): x1, f(x).diff(x).subs(x, x2): x3}
ans = dsolve(ode,y, ics={y.subs(x,0):1,y.diff(x).subs(x,0):2})
print('二階ode的解：', ans.args[0],'=',ans.args[1])

#畫出所解出的原函數解
#plot(ans.args[1],line_color='red')
plot(ans.args[1],(x,-1,1),line_color='red')

成果圖示：
成果圖示：
成果圖示：
第一種微分寫法程式碼
第二種微分寫法程式碼

範例8-2：二階常係數齊性ODE的解法(二)--重根

2.範例8-2：二階常係數齊性ODE的解法(二)--重根
# ode ： y''(x) + 4y'(x) +4y = 0
#邊界條件或起始條件：y(0) =1, y'(0) =2
#語法：dsolve(eq, f(x), ics={y(0):1,y(x).diff(x).subs(x,0):2})
#所解出的原方程式： y(x) = (4*x + 1)*exp(-2*x)
from sympy import *
x = symbols('x')
y = Function('y')(x)
#注意：ode是個方程式，不是函數
ode = Eq(y.diff(x,2) + 4*y.diff(x) + 4*y, 0)
ans = dsolve(ode, y, ics={y.subs(x,0):1, y.diff(x).subs(x,0):2})
print('二階ode的解：', ans.args[0],'=',ans.args[1])

#畫圖y(x)
plot(ans.args[1],(x,-1, 5), line_color='red')

成果圖示：
成果圖示：
第一種微分寫法程式碼
第二種微分寫法程式碼

3.範例8-3：二階常係數齊性ODE的解法(三)--複數根

3.範例8-3：二階常係數齊性ODE的解法(三)--複數根
# ode ： y''(x) + 2y'(x) + 5y = 0
#邊界條件或起始條件：y(0) =1, y'(0) =1
#語法：dsolve(eq, f(x), ics={y(0):1,y(x).diff(x).subs(x,0):2})
#所解出的原方程式：y(x) = (sin(2*x) + cos(2*x))*exp(-x)
from sympy import *
x = symbols('x')
y = Function('y')
ode = diff(y(x),x,2) + 2*diff(y(x),x) + 5*y(x)
ans = dsolve(ode, y(x), ics={y(0):1, y(x).diff(x).subs(x,0):1})
print('二階ode的解：',ans.args[0],'=',ans.args[1])

成果圖示：
成果圖示：
第一種微分寫法程式碼
第二種微分寫法程式碼

範例8-4：Euler-Cauchy 方程式之解--複數根

4.範例8-4：Euler-Cauchy 方程式之解--複數根
# ode ：x**2*y''(x) - x*y'(x) + 4y = 0
#邊界條件或起始條件：y(1) =2, y'(1) =-1
#語法：dsolve(eq, f(x), ics={y(1):2,y(x).diff(x).subs(x,1):-1})
#所解出的原方程式：y(x) = -sin(2*log(x))/2 + 2*cos(2*log(x))
from sympy import *
x = symbols('x')
y = Function('y')(x)
#注意：ode是個方程式，不是函數
ode = Eq(x**2*y.diff(x,2) + x*y.diff(x) + 4*y)
ans = dsolve(ode, y, ics={y.subs(x,1):2, y.diff(x).subs(x,1):-1})
print('ode的解：', ans.args[0], '=', ans.args[1])

plot(ans.args[1], (x,0,10), line_color='red')

成果圖示：
成果圖示：
第一種微分寫法程式碼
第二種微分寫法程式碼

範例8-5：二階常係數非齊性ODE之特解

5-範例8-5：二階常係數非齊性ODE之特解：y''(x)+2y'(x)-3y(x)=3x**2，邊界條件為y(0)=0, y'(0)=0
#因為多一個4x**2，這個是非齊性項（全部的係數都是固定常數）
#ode ：y''(x)+2y'(x)-3y(x)=3x**2
#邊界條件或起始條件：y(0)=0, y'(0)=0
#語法：dsolve(eq, f(x), ics={y(0):0,y(x).diff(x).subs(x,0):-0})
#所解出的原方程式：y(x) = -x**2 - 4*x/3 + 3*exp(x)/2 - 14/9 + exp(-3*x)/18
from sympy.interactive import printing
printing.init_printing(use_latex=True)
from sympy import *
x = symbols('x')
y = Function('y')(x)
#注意：ode是個方程式，不是函數
ode = y.diff(x,2)+2*y.diff(x)-3*y - 3*x**2
ans = dsolve(ode, y, ics={y.subs(x,0):0, y.diff(x).subs(x,0):0})
print('二階常係數非齊性ODE之特解：',ans.args[0],'=',ans.args[1])

#畫圖y(x)
plot(ans.args[1],(x,-10,0), line_color='red')

成果圖示：
成果圖示：
第一種微分寫法程式碼
第二種微分寫法程式碼

範例8-6：高階(三階)常係數齊性 ODE 之求解

6.範例8-6：高階(三階)常係數齊性 ODE 之求解：y''''(x)-y''(x)-y'(x)+y(x)=0，邊界條件為y(0)=2, y'(0)=1, y''(0)=0
#因為多一個4x**2，這個是非齊性項（全部的係數都是固定常數）
#ode ：y''''(x)-y''(x)-y'(x)+y(x)=0
#邊界條件或起始條件：y(0)=2, y'(0)=1, y''(0)=0
#語法：dsolve(eq, f(x), ics={y(0):0,y(x).diff(x).subs(x,0):-0})
#所解出的原方程式：y(x) = -x**2 - 4*x/3 + 3*exp(x)/2 - 14/9 + exp(-3*x)/18
from sympy import *
x = symbols('x')
y = Function('y')
#注意：ode是個方程式，不是函數
ode = Eq(diff(y(x),x,3)-diff(y(x),x,2)-diff(y(x),x)+y(x),0)
ans = dsolve(ode,y(x),ics={y(0):2, y(x).diff(x).subs(x,0):1, y(x).diff(x,2).subs(x,0):0})
#ans = dsolve(ode,y(x),ics={y(0):2},y(x).diff(x).subs(x,0):1, y(x).diff(x).diff(x).subs(x,0):0)
print('高階常係數齊性 ODE 之求解：',ans.args[0],'=',ans.args[1])

#顯示數學符號函數
from sympy.interactive import printing
printing.init_printing(use_latex=True)
display(y)

#畫圖y(x)
plot(ans.args[1],(x,-10,10), line_color='red')

成果圖示：
成果圖示：
第一種微分寫法程式碼
第二種微分寫法程式碼

範例8-7：一階聯立微分方程式，求特解

7.範例8-7：一階聯立微分方程式，求特解
#eq1：y1'(t) + 3y1 = 5
#eq1：y2'(t) - 3y2 = 5
#起始條件：y1(0)=1, y2(0)=2
#關鍵方法：dsolve()本身就可以解聯立微分方程式組
#語法：dsolve(聯立方程式list, [y1,y2]list)
#邊界條件：ics={y1.subs(t,0):1, y2.subs(t,0):2}
from sympy import *
t = symbols('t')
y1 = Function('y1')(t)
y2 = Function('y2')(t)
ode1 = Eq(y1.diff(t) + 3*y1, 5)
ode2 = Eq(y2.diff(t) - 3*y2, 5)
ans = dsolve([ode1, ode2], [y1, y2], ics={y1.subs(t,0):1, y2.subs(t,0):2})
print('y1=', ans[0].rhs)
print('y2=', ans[1].rhs)

#畫圖
f1 =ans[0].rhs
f2 =ans[1].rhs
#print(f)
p1 = plot(f1, f2, (t,-1,0.5), show=False)
#設定顏色
p1[0].line_color = 'r'
p1[1].line_color = 'g'
#設定圖例legend
p1[0].label = '$y1$'
p1[1].label = '$y2$'
#要打開 True
p1.legend = True
p1.show()

成果圖示：
第一種微分寫法程式碼
第二種微分寫法程式碼

範例8-8：二階聯立微分方程式，求通解(y1,y2)，聯立齊性 ODE 的解法

8.範例8-8：二階聯立微分方程式，求通解(y1,y2)，聯立齊性 ODE 的解法
#eq1：y1''(x) = -5y1 + 2y2
#eq1：y2''(x) = 2y1 - 2y2
#關鍵方法：dsolve()本身就可以解聯立微分方程式組
#語法：dsolve(聯立方程式list, [y1,y2]list)
#ans = dsolve([[ode1, ode2]], [y1,y2])
from sympy import *
x = symbols('x')
y1 = Function('y1')(x)
y2 = Function('y2')(x)
ode1 = Eq(y1.diff(x,2), -5*y1 + 2*y2)
ode2 = Eq(y2.diff(x,2), 2*y1 - 2*y2)
#解聯立方程式的方法：先把所有的方程式放入一個eq 陣列
eqlis = [ode1, ode2]
#dsolve(聯立方程式list, [y1,y2]list)
ans = dsolve(eqlis, [y1,y2])
#print(ans)

#顯示標準數學符號
from sympy.interactive import printing
printing.init_printing(use_latex=True)
display(ans)
#ans有兩組解答，如何取得：第一個=ans[0], 第二個=ans[1]
display(ans[0])
display(ans[1])
# expr = Eq（方程式，equation) 如何取得其中的函數部分
#expr.lhs：取得方程式的左邊
#expr.rhs：取得方程式的右邊
print('y1=',ans[0].rhs)
print('y2=',ans[1].rhs)

成果圖示：
第一種微分寫法程式碼

1.聯立方程式的矩陣求法

1.聯立方程式的矩陣求法

範例9-1：聯立方程式的聯立求法（線性代數）

2.範例9-1：聯立方程式的聯立求法（線性代數）
from sympy import *
from sympy.interactive import printing
printing.init_printing(use_latex=True)

x, y = symbols('x y')
#二元一次聯立方程式
eq1 = Eq(x + 2*y - 8, 0)
eq2 = Eq(2*x - y - 6, 0)
row1 = [1,2,8]
row2 = [2,-1,6]
#注意：格式：Matrix((row1, row2))
system = Matrix((row1, row2))
display(system)

#1.第一種：使用線性代數的函數，解聯立方程式
print('1.第一種：使用線性代數的函數，解聯立方程式:')
#格式：solve_linear_system(矩陣,變數x,變數y...)
ans = solve_linear_system(system,x,y)
print(ans)
#注意：這個資料格式{x: -4, y: -2}，像dict，但不是，因為若是dict會是{'x': -4, 'y': -2}
#dict的取值：ans['x']
#目前資料格式的取法ans[x]
print('x=',ans[x])
print('y=',ans[y])

#2.第二種：直接用solve就可以解聯立方程式
print('2.第二種：直接用solve就可以解聯立方程式:')
ans = solve([eq1,eq2],[x,y])
print('x=',ans[x])
print('y=',ans[y])

#畫圖
p1 = None
p1 = plot_implicit(eq1,line_color='red',how=False)
p2 = plot_implicit(eq2,line_color='blue',show=False)
p1.extend(p2)
p1.show()

成果圖示：
程式碼內容

1.泰勒展開：以多項式來逼近函數

1.泰勒展開：以多項式來逼近函數
泰勒級數是由牛頓的學生Taylor所建立的定理
（1）特色：任何的多項式，都可以寫成以a為參考點的多項式（x-a)^n+(x-a)^n-1+.....

（2）例如1：x**3 + x**2 + x +1，可以轉換成以1的參考點（x-1）的多項式
利用綜合除法可以得到
x**3 + x**2 + x +1
= 4 + 6(x-1) + 4(x-1)**2 + (x-1)**3
=Taylor Series習慣用升冪級數

（3）例如1：2x**3 - x**2 - x* +1，可以轉換成以-1的參考點（x+1）的多項式
利用綜合除法可以得到
2x**3 - x**2 - x* +1
= -1 + 7(x+1) - 7(x+1)**2 + 2(x+1)**3

（4）一個n次方的多項式


f(x) = An*x**n + An-1*x**(n-1) + An-2*x**(n-2)+...+ A0
若以b為參考點（x-b），可以轉換成泰勒級數：
=C0 + C1(x-b) + C2(x-b)**2 + C3(x-b)**3 +....Cn(x-b)**n
其中：
C0 = f(b)
Cn = An

（5）泰勒形式（通式）：以a為參考點的泰勒形式
= =C0 + C1(x-a) + C2(x-a)**2 + C3(x-a)**3 +....Cn(x-a)**n

（5）泰勒形式（標準式）：以 0 為參考點的泰勒形式
= =C0 + C1x + C2x**2 + C3x**3 +....Cnx**n

2.泰勒級數的應用-1：大約估值

2.泰勒級數的應用-1：大約估值
（1）例如：某函數f(x) = x**3 + x**2 + x + 1，請估計當f(0.02)的估計值（誤差到小數點第二位）
轉成以0為參考點的泰勒級數
f(x) = C0 + C1(x-0) + C2(x-0)**2 + ....Cn(x-0)**n
= 1 + x + x**2 + x**3
= 1 + 0.02 + 0.02**2 + 0.02**3
因為只要到小數點第二位，所以0.0**2以後的誤差都可以不用考慮
= 1 + 0.02
估計值 = 1.02

（2）例如：f(x) = x**3 - x**2 - x -1，請估計f(1.97)的大概值到小數點第二位
f(x) = x**3 - x**2 - x -1
轉成(x-2)的泰勒級數
= C1 + C2(x-2) +不考慮後面的誤差
= 1 + 7(x-2)
= 1 + 7 * -0.03 = 1-0.21 = 0.79

因為其x=2的估計值 f(2) 接近= 1+7（x-2) = 固定斜率的直線
所以
若是畫出f(x)的圖形
在x=2出的圖形附近區間，f(x)是非常逼近直線的

3.泰勒一次項係數即導數

3.泰勒一次項係數即導數
（1）泰勒形式（通式）：以a為參考點的泰勒形式
= =C0 + C1(x-a) + C2(x-a)**2 + C3(x-a)**3 +....Cn(x-a)**n
其中誤差到小數點第1-2位的是在一次項係數 = C1

（2）一次項係數C1就是f(a)的導數
C1 = f'(a) = f在a的導數

4.用微分計算泰勒多項式的係數

4.用微分計算泰勒多項式的係數
（1）泰勒形式（通式）：以a為參考點的泰勒形式
= C0 + C1(x-a) + C2(x-a)**2 + C3(x-a)**3 +....Cn(x-a)**n
= f(a)+
f'(a)(x-1)+
f''(a)/2*(x-a)**2 +
f'''(a)/6*(x-a)**3 +
..
f''''''(a)/n!*(x-a)**n

其中的係數規律：
C0 = f(a)
C1 = f'(a)
C2 = 1/2*f''(a)
Cn = An
Ck = (1/k!)*f''''(a)，k次導數

（2）應用1：f(x) = x**3 + x**2 + x + 1
求f(0.97)的估計值，與精確值

A).估計值：
f(x) = x**3 + x**2 + x + 1
轉成泰勒級數
= C0 + C1(x-1) .....不用考慮誤差
= f(a) + f'(a)(x-1)
= f(1) + f'(1)(x-1)
= 4 + 6(0.03) = 4.12
其中 f'(1) = 3x**2 + 2x +1 = 3+2+1 = 6

B).精確值：
f(x) = x**3 + x**2 + x + 1
轉成泰勒級數
= C0 + C1(x-1) ....C3(x-1)**3
= f(a) + f'(a)(x-1) + 1/2!*f''(a)(x-1)**2 + 1/3!*f'''(a)(x-1)**3
= f(1) + 6(x-1) + 8/2*(x-1)**2 + 1/3*(x-1)**3
= 4 + 6(-0.03) + 4(-0.03)**2 + 1*-0.03**3
= 4 -0.18 +0.0036 - 0.000027
= 3.823573
其中
C1 = f'(1) = 3x**2 + 2x +1 = 3+2+1 = 6
C2 = 1/2*f''(1) = 1/2(6x + 2) = 4
C 3 = f'''(1) = 1

5.多項式函數的圖形特色

5.多項式函數的圖形特色
（1）多項式函數的特徵
二次多項式：必定有極大值，極小值（但不一定有反曲點）
三次多項式：必定有反曲點（不一定有極值）
n次多項式（n>3)：可能有反曲點，若有反曲點，必須要f''(x)二階導數=0

（2）三次多項式圖形特徵：必定有反曲點

（3）n次多項式（n>3)：若二階導數f''(x)=0，就會有反曲點

（4）反曲點的應用1：例如股票線型，若是能夠發現出現反曲點了，表示曲線將出現凹口向上，或凹口向下

（5）如何求得三次多項式圖形的局部極大值，極小值
A).泰勒級數的一次導數係數f'(x) = C1 = 0，此時的x0就是極值（有兩個解，分別是極大值，與極小值）
B）.如何判別極大值，或是極小值
若二階導數f''(x0) >0，為正，那麼表示曲線開口向上，代表，x0處是極小值
若二階導數f''(x0) <0，為負，那麼表示曲線開口向下，代表，x0處是極大值

6.典型範例：由銷售量q與單價x的三次多項式關係，求月營收的極大值極小值

6.典型範例：由銷售量q與單價x的三次多項式關係，求月營收的極大值極小值，分別是在什麼的定價時發生的？
（1）銷售量q與訂價x的關係
q = 1/5*(x-290)(x-300)-25(x-300)+200
月營收f = 銷售量q *定價x = qx
f = qx = x*(1/5*(x-290)(x-300)-25(x-300)+200)
f = 1/5*x**3-143x**2+26900x = 三次多項式

（2）若有極值發生，是在一階導數係數C1=0時發生的兩個解
f'(x) = C1 = 0
3/5**2 - 286x +26900 = 0
兩個解
x1 = 715/3 + 5*sqrt(4309)/3 = 128.9..
x2 = 715/3 - 5*sqrt(4309)/3 = 347.7..

（3）判別是極大值，或是極小值？
看反曲點位置，判別開口向上，或開口向下？
f''(x) = 0
6/5* x -286 =0
A).在123.9位置
f''(128.9..) = 6/5*128.9 -286 < 0
開口向下，
在x1 = 128.9處，發生極大值

A).在347.7位置
f''(347.7..) = > 0
開口向上，
在x1 = 347.7處，發生極小值

7.完整的n次泰勒多項式函數（以導數表示係數），及其發生極值的位置

7.完整的n次泰勒多項式函數（以導數表示係數），及其發生極值的位置
（1）n次多項式函數的泰勒級數（以導數表示係數）





（2）多項式判別判別極值的方法：
#A).泰勒級數的一次導數係數f'(x-x0) = C1 = 0，此時的x0就是極值（有兩個解，分別是極大值，與極小值）
#B）.如何判別極大值，或是極小值
#若二階導數f''(x0) >0，為正，那麼表示曲線開口向上，代表，x0處是極小值
#若二階導數f''(x0) <0，為負，那麼表示曲線開口向下，代表，x0處是極大值
from sympy import *

範例10-1：使用sympy把多項式轉成泰勒級數

8.範例10-1：使用sympy把多項式轉成泰勒級數
#多項式函數f(x) = x**3 + x**2 + x* +1
#以1為參考點（x-1），轉成泰勒級數
from sympy import *
from sympy.interactive import printing
printing.init_printing(use_latex=True)

x = symbols('x')
f = Function('f')(x)
f = x**3 + x**2 + x +1
#轉換成級數series
#(1).第一種語法：expr.series(x=None, x0=0, n=6, dir='+')
#(2).第二種語法：series(expr, x=None, x0=0, n=6, dir='+')
#缺點，會產生O項目，還要就經過.removeO()，較麻煩

#(1).第一種語法：expr.series(x=None, x0=0, n=6, dir='+')
tysr = f.series( x, x0=1, n=3)
display(tysr)
#把O項刪除
tysr = f.series(x, x0=1, n=3).removeO()
display(tysr)
#最高次方的係數為O, 故修改n3 = 4
tysr =f.series(x, x0=1, n=4).removeO()
display(tysr)

#(2).第二種語法：series(expr, x=None, x0=0, n=6, dir='+')
tysr =series(f, x, x0=1, n=4)
display(tysr)

#缺點：但是6x沒有轉換成(x-1)的格式

#(3).f(0.02)的估計值1.02，請算精確值
valueoftysr = tysr.subs(x,0.02)
valueoff = f.subs(x, 0.02)
print('泰勒級數在x=0.02的計算結果=',valueoftysr)
print('f(x)在x=0.02的計算結果=',valueoff)
print('用泰勒級數一次項的大概估計值=',1.02)

#畫出兩條圖
p1 = plot(f,tysr,(x,-10,10))
p1[0].line_color = 'red'
p1[1].line_color = 'blue'
p1[0].label = '$f(x)$'
p1[1].label = '$taylor series$'
p1.legend = True
p1.show()

成果圖示：
程式碼內容

範例10-2：f(1.97)泰勒級數估計值大概為-76.03，請計算精確值

9.範例10-2：f(1.97)泰勒級數估計值大概為-76.03，請計算精確值
#多項式函數: f(x) = x**3 - x**2 - 40x -1
#以2為參考點（x-2），轉成泰勒級數
from sympy import *
x = symbols('x')
f = Function('y')(x)
f = x**3 -x**2 -40*x -1

ty1 = series(f, x, x0=2, n=4)
print('泰勒級數的精確值=', ty1.subs(x,1.97))
print('f(x=2)的精確值=', f.subs(x,1.97))
print('泰勒級數用一次項的估計值=', -76.03)

#印出標準方程式
from sympy.interactive import printing
printing.init_printing(use_latex=True)
display(f)
display(ty1)

#畫出兩條圖
p1 = plot(f,ty1)
p1[0].line_color = 'red'
p1[1].line_color = 'blue'
p1[0].label = '$f(x)$'
p1[1].label = '$taylor series$'
p1.legend = True
p1.show()

成果圖示：
程式碼內容

範例10-3：求三次多項式的極值發生地方

10.範例10-3：求三次多項式的極值發生地方
#多項式函數: f(x) = x**3 - x**2 - x -1
#多項式判別判別極值的方法：
#A).泰勒級數的一次導數係數f'(x-x0) = C1 = 0，此時的x0就是極值（有兩個解，分別是極大值，與極小值）
#B）.如何判別極大值，或是極小值
#若二階導數f''(x0) >0，為正，那麼表示曲線開口向上，代表，x0處是極小值
#若二階導數f''(x0) <0，為負，那麼表示曲線開口向下，代表，x0處是極大值
from sympy import *
from sympy.interactive import printing
printing.init_printing(use_latext=True)
x = symbols('x')
f = Function('f')(x)
f = x**3 - x**2 - x -1
#先算f'(x)一階導數
f1 = f.diff(x)
eqf1 = Eq(f1, 0)
ans = solve(eqf1, x)
x1 = ans[0]
x2 = ans[1]
print('解x1=', x1)
print('解x2=', x2)
#由反曲點=f''(x)判別開口向上，向下？
f2 = f.diff(x,2)
if f2.subs(x, x1) > 0:
print('x1=',x1,'處，是極小值')
else:
print('x1=',x1,'處，是極大值')

if f2.subs(x, x2) > 0:
print('x2=',x2,'處，是極小值')
else:
print('x2=',x2,'處，是極大值')

#畫出圖形驗證
plot(f,(x, -2, 2), line_color='red')

成果圖示：
程式碼內容

例10-4：由銷售量q與單價x的三次多項式關係，求月營收的極大值極小值

11.範例10-4：由銷售量q與單價x的三次多項式關係，求月營收的極大值極小值，分別是在什麼的定價時發生的？
#（1）月營收f 與產品定價的關係式：三次多項式
#f(x) = 1/5*x**3-143x**2+26900x
#多項式判別判別極值的方法：
#A).泰勒級數的一次導數係數f'(x-x0) = C1 = 0，此時的x0就是極值（有兩個解，分別是極大值，與極小值）
#B）.如何判別極大值，或是極小值
#若二階導數f''(x0) >0，為正，那麼表示曲線開口向上，代表，x0處是極小值
#若二階導數f''(x0) <0，為負，那麼表示曲線開口向下，代表，x0處是極大值
from sympy import *
from sympy.interactive import printing
printing.init_printing(use_latex=True)

x = symbols('x')
f = Function('f')(x)
f = (1/5)*x**3 - 143*x**2 + 26900*x
#判別月營收f發生極值的dif：f'(x-x0) =0
f1 = f.diff(x)
eqf1 = Eq(f1, 0)
ans = solve(eqf1, x)
x1 = ans[0]
x2 = ans[1]
#判別極值處，是極大值，，或是極小值（看f''(x-x0)正負，判別開口向上或向下
f2 = f.diff(x,2)
if f2.subs(x, x1) > 0:
print('x1=',x1,'處，有極小值')
else:
print('x1=',x1,'處，有極大值')
if f2.subs(x, x2) > 0:
print('x2=',x2,'處，有極小值')
else:
print('x2=',x2,'處，有極大值')

#畫圖驗證
plot(f, (x, 0, 450), line_color='red')

成果圖示：
程式碼內容

2.正餘弦sin(x),cos(x)的泰勒多項式

12.正餘弦sin(x),cos(x)的泰勒多項式
（1）觀念：任何的函數都可以用泰勒多項式來逼近
理論上只要n次方取的夠多，就能夠逼近原函數
（但是在實際上，有些函數圖形，無法以泰勒級數來逼近，例如：對數函數ln(x)）

（2）以泰勒級數（參考點x0=0）來表示 sin(x)，
重點：只要知道泰勒多項式的各個係數Cn，那麼就可以逼近函數sin(x)了
f'(x) = sin(x)
f''(x) = cos(x)
f'''(x) = -sin(x)
f''''(x) = -cos(x)
f'''''(x) = sin(x)
f''''''(x) = cos(x)
若是x0=0，則
f'(0) = 0
f''(0) = 1
f'''(0) = 0
f''''(0) = -1
f'''''(0) = 0
f''''''(0) = 1

Sn = x -1/3!*x**3 + 1/5!*x**5 -1/7!*x**7 +1/9!*x**9 -1/11!*x**11 +....
（偶次方的係數都是0，只剩下奇次方)

當n次方越大，就會越來越接近sin(x)

（3）以泰勒級數（參考點x0=0）來表示 cos(x)，
因為sin'(x) = cos(x)
而
Sn = x -1/3!*x**3 + 1/5!*x**5 -1/7!*x**7 +1/9!*x**9 -1/11!*x**11 +....
Cn = [Sn+1]' =
Cn = 1 - 1/2!*x**2 +1/4!*x**4 -1/6!*x**6 +1/8！*x**8 -1/10!*x**10 +..
只剩下偶次方)

13.指數函數的泰勒多項式

13.指數函數的泰勒多項式
（2）以泰勒級數（參考點x0=0）來表示標準指數函數 exp(x)
重點：只要知道泰勒多項式的各個係數Cn，那麼就可以逼近函數exp(x)了
f'(x) = exp(x)
f''(x) = exp(x)
f'''(x) = exp(x)
f''''(x) = exp(x)
f'''''(x) = exp(x)
f''''''(x) = exp(x)
若是x0=0，則exp(0) = 1
f'(0) = 1
f''(0) = 1
f'''(0) = 1
f''''(0) = 1
f'''''(0) = 1
f''''''(0) = 1

結論：e**x若以0為參考點，的n次泰勒多項式為
En = 1 + x 1/2!*x**2 + 1/3!*x**3 +1/4!*x**4 +1/5!*x**5 +1/6!*x**6 +....
（偶次方的係數都是0，只剩下奇次方)

當n次方越大，就會越來越接近exp(x)

13.對數函數的泰勒多項式

13.對數函數的泰勒多項式
（1）以泰勒級數（參考點x0=1）來表示標準對數函數 ln(x)
因為0不在ln(x)的定義範圍內，所以用x0=1
重點：只要知道泰勒多項式的各個係數Cn，那麼就可以逼近函數ln(x)了
若以1為參考點，的n次泰勒多項式為
En = (x-1) - 1/2*(x-1)**2 + 1/3*(x-1)**3 -1/4*(x-1)**4
（偶次方的係數都是0，只剩下奇次方)

（2）但是對數的泰勒級數逼近原函數的效果較差，比上述的sin(x),exp(x)都發生較大的誤差
但是在x0 = a = 2的地方，是逼近原函數的，但是在遠離x0=2的地方，誤差就較大了

（3）在實際上，有些函數圖形，無法以泰勒級數來逼近，例如：對數函數ln(x)

範例10-5：以泰勒多項式逼近指數函數exp(x)

14.範例10-5：以泰勒多項式逼近指數函數exp(x)
from sympy import *
x = symbols('x')
f = Function('f')(x)
f = exp(x)
#畫圖比較
ty2 = series(f, x, x0=0, n=3).removeO()
ty4 = series(f, x, x0=0, n=5).removeO()
ty8 = series(f, x, x0=0, n=9).removeO()

p1 = plot(f, ty2, ty4, ty8, (x, -4, 4), show=False)
p1[0].line_color = 'red'
p1[1].line_color = 'blue'
p1[2].line_color = 'green'
p1[3].line_color = 'yellow'
p1[0].label = '$e(x)$'
p1[1].label = '$taylor(n=2)$'
p1[2].label = '$taylor(n=4)$'
p1[3].label = '$taylor(n=8)$'
p1.legend = True
p1.show()

成果圖示：
程式碼內容

結果：使用n=8次方泰勒多項式時，就非常逼近exp(x)了

範例10-6：以泰勒多項式逼近sin(x)函數

15.範例10-6：以泰勒多項式逼近sin(x)函數
from sympy import *
x = symbols('x')
f = Function('f')(x)
f = sin(x)

ty1 = series(f, x, x0=0, n=24).removeO()
ty2 = series(f, x, x0=0, n=25).removeO()
ty3 = series(f, x, x0=0, n=26).removeO()
#畫圖
p1 = plot(f, ty1, ty2, ty3, (x, -10, 10), show='False')
p1[0].line_color= 'red'
p1[1].line_color= 'blue'
p1[2].line_color= 'green'
p1[3].line_color= 'yellow'
p1[0].label= '$sin(x)$'
p1[1].label= '$taylor(n=23)$'
p1[2].label= '$taylor(n=24)$'
p1[3].label= '$taylor(n=25)$'
p1.legend = True
p1.show()

成果圖示：
程式碼內容
結果：使用n=25次方泰勒多項式時，就非常逼近sin(x)了 (x在-10~10的範圍內）

範例10-7：以泰勒多項式逼近ln(x)函數

16.範例10-7：以泰勒多項式逼近ln(x)函數
#（1）以泰勒級數（參考點x0=1）來表示自然對數函數 ln(x)
#因為0不在ln(x)的定義範圍內，所以用x0=1
#En = (x-1) - 1/2*(x-1)**2 + 1/3*(x-1)**3 -1/4*(x-1)**4
from sympy import *
x = symbols('x')
f = Function('f')(x)
# 在sympy裡面，自然對數指令：ln(x) = log(x)
f = ln(x)

ty1 = series(f, x, x0=1, n=2).removeO()
ty2 = series(f, x, x0=1, n=3).removeO()
ty3 = series(f, x, x0=1, n=5).removeO()
#若是設定n=30，則圖形更為失真，造成全部的圖形都無法顯示出來
#總之，泰勒級數無法逼近ln(x)的指數圖形
#畫圖
p1 = plot(f, ty1, ty2, ty3, (x, -1, 2), show='False')
p1[0].line_color= 'red'
p1[1].line_color= 'blue'
p1[2].line_color= 'green'
p1[3].line_color= 'yellow'
p1[0].label= '$sin(x)$'
p1[1].label= '$taylor(n=23)$'
p1[2].label= '$taylor(n=24)$'
p1[3].label= '$taylor(n=25)$'
p1.legend = True
p1.show()

成果圖示：
程式碼內容
結論：
（1）以泰勒級數逼近ln(x)的指數圖形，最後失敗
（2）所以不是所有的函數圖形，都難以泰勒級數來逼近

範例14-9：讀取網址的網頁

1.傅立葉級數的應用

1.傅立葉級數的應用
在工程上使用傅立葉級數的機會很多。
例如：訊號處理：電磁波變化，溫度訊號，心電圖，腦波
音響，音效波形，聲音辨識

2.傅立葉數學的學習內容與步驟

2.傅立葉數學的學習內容與步驟
（1）傅立葉級數（有兩種：週期為為2𝝿，或是2L)
（2）傅立葉積分（週期為∞)，傅立葉cosine積分，sine積分
（3）傅立葉轉換（週期為∞),
（4）Laplace轉換（需要調整積分函數型態）

3.傅立葉級數公式

3.傅立葉級數公式
（1）一個波函數f(x)，可以用傅立葉級數公式來模擬之：
（2）傅立葉級數有兩種寫法，
一種是以sin(x)+cos(x)合成波
一種是以諧波exp(inπx)的複數指數函數的和合成波


（3）傅立葉級數, sine(2πnx)-cos(2πnx) 形式：
週期為 2π的函數 f(x) 之 Fourier 級數



（4）傅立葉級數, 以諧波exp(inπx)的複數指數函數的和合成波形式：


（5）原始函數f(t)的相位與頻率
如果原始波是：sin(θ)
不同頻率波：sin(θ/f)，例如：sin(θ)，sin(2θ)，sin(3θ)，sin(4θ)...就是不同頻率f的波
不同相位波：sin(θ+φ)，sin(θ+2φ)+sin(θ+3φ)....就是不同相位φ的波

4.傅立葉級數的物理意義

4.傅立葉級數的物理意義：
（1）由傅立葉級數（以諧波exp(inπx)的複數指數函數的和合成波）形式來看：

代表意義：任何一個週期波f(x)，都可以由很多不同頻率nπ的諧波（單頻）組合而成。

（2）什麼是諧波？
諧波：就是會產生共振的波。

什麼波才能共振？
必須要是整數倍頻率的波，才能共振。
例如：exp(inπx) + exp(i2nπx) + exp(i3nπx) + exp(i4nπx)

（3）exp(inπx)是世界上唯一的單頻率的波
所以某個波函數f(x)，就可以由這些單頻波exp(inπx)組合而成。

（3）範例：用傅立葉級數模擬鋸齒波
我們現在用上面的公式給出一個簡單函數的傅立葉級數展開式。考慮一個鋸齒波

我們用5個不同倍率的諧波組合成傅立葉級數來模擬鋸齒波
exp(inπx) + exp(i2nπx) + exp(i3nπx) + exp(i4nπx)+ exp(i5nπx)




其中的Cn：代表第n個諧波exp(inπx) 的分量
其中的C0：代表週期內單位時間的平均值，代表會把整個波上移或下移的平均值。
（在電壓波，也是直流電的電壓值）

（4）範例：用傅立葉級數（n=5個不同頻率但成倍數的諧波exp(inπx)來模擬方波


（5）角頻率ω，與弧頻率ξ，時間週期T的關係
ω = 2πξ
ω = 2π/T
（T為時間週期）

5.傅立葉級數公式的視覺化Series

5.傅立葉級數公式的視覺化
（1）傅立葉級數公式的視覺化
https://www.youtube.com/watch?v=r18Gi8lSkfM
https://www.youtube.com/watch?v=spUNpyF58BY
https://www.youtube.com/watch?v=spUNpyF58BY&t=938s

（2）可以用傅立葉級數模擬複雜圖形，例如：人的大腦圖形
影片：
https://www.youtube.com/watch?v=ds0cmAV-Yek

（3）熱傳非線性偏微分方程式，也可用傅立葉級數來模擬原始函數T溫度變化的過程
影片：
https://www.youtube.com/watch?v=r6sGWTCMz2k
https://www.youtube.com/watch?v=ToIXSwZ1pJU

範例11-1：在(-pi < x < pi)內用傅立葉級數模擬三次方多項式的線型

6.範例11-1：在(-pi < x < pi)內用傅立葉級數Fourier Series模擬三次方多項式的線型
from sympy import *
from sympy.interactive import printing
printing.init_printing(use_latex=True)

x = symbols('x')
f = Function('f')(x)
#f = x**3 + x**2 + x +1
f = x**3 -x**2 -40*x -1
#語法：fourier_series(要模擬的f函數, (x的範圍, -pi, pi))
s = fourier_series(f, (x, -pi, pi))
#不太可能去使用n=∞去計算，在實際計算時要截短級數truncate（指定n=數字）
fs5 = s.truncate(n=5)
display(fs5)
plot(fs5)

#再模擬另外一種線型的三次方函數
f = x**3 -x**2 -40*x -1
s = fourier_series(f, (x, -pi, pi))
fs5 = s.truncate(n=5)
plot(fs5)

成果圖示：
程式碼內容

範例11-2：比較不同n級數的傅立葉級數去模擬三次方多項式的

7.範例11-2：比較不同n級數的傅立葉級數Fourier Series去模擬三次方多項式的線型效果
#（1）若是n級數太少(n=5)，就會形成鋸齒狀模擬結果
#（2）若是n級數夠大(n=40))，就會形成平滑化的模擬結果
from sympy import *
from sympy.interactive import printing
printing.init_printing(use_latex=True)

x = symbols('x')
f = Function('f')(x)
f = x**3 + x**2 + x +1
#語法：fourier_series(要模擬的f函數, (x的範圍, -pi, pi))
s = fourier_series(f, (x, -pi, pi))
#不太可能去使用n=∞去計算，在實際計算時要截短級數truncate（指定n=數字）
fs5 = s.truncate(n=5)
fs10 = s.truncate(n=40)

#畫圖
p1 = plot(fs5,fs10,show=False)
p1[0].line_color = 'red'
p1[1].line_color = 'green'
p1[0].label = '$fourier(n=5)$'
p1[1].label = '$fourier(n=40)$'
p1.legend = True
p1.show()

成果圖示：
程式碼內容

範例11-3：比較不同n級數的傅立葉級數去模擬指數函數exp(x)

8.範例11-3：比較不同n級數的傅立葉級數Fourier Series去模擬指數函數exp(x)的線型效果
#（1）若是n級數太少(n=5)，就會形成鋸齒狀模擬結果
#（2）若是n級數夠大(n=50))，就會形成平滑化的模擬結果
from sympy import *
x = symbols('x')
f = Function('f')(x)
f = exp(x)
s = fourier_series(f, (x, -pi, pi))
fs5 = s.truncate(n=5)
fs50 = s.truncate(n=50)

p1 = plot(fs5, fs50, show=False)
p1[0].line_color = 'red'
p1[1].line_color = 'green'
p1[0].label = '$fourier(n=5)$'
p1[1].label = '$fourier(n=50)$'
p1.legend = True
p1.show()

成果圖示：
程式碼內容

範例11-4：比較不同n級數的傅立葉級數去模擬自然對數函數ln(x)

9.範例11-4：比較不同n級數的傅立葉級數Fourier Series去模擬自然對數函數ln(x)的線型效果
#注意：若是要傅立葉級數Fourier Series去模擬自然對數函數ln(x)的線型，因為，自然對數函數ln(x)在0的位置的值為-∞（負無限大），所以在x=0位置無法模擬
#本題採用 (0.5pi < x < 2.5pi)的週期波
#（1）若是n級數太少(n=5)，就會形成鋸齒狀模擬結果
#（2）若是n級數夠大(n=40))，就會形成平滑化的模擬結果（但是計算時間要等久一點）
from sympy import *
x = symbols('x')
f = Function('f')(x)
f = ln(x)
s = fourier_series(f, (x, pi/2, 5/2*pi))
fs5 = s.truncate(n=5)
fs50 = s.truncate(n=40)

p1 = plot(fs5, fs50, show=False)
p1[0].line_color = 'red'
p1[1].line_color = 'green'
p1[0].label = '$fourier(n=5)$'
p1[1].label = '$fourier(n=40)$'
p1.legend = True
p1.show()

成果圖示：
程式碼內容

範例11-5：求週期函數Square Wave方波函數f(x)，在(-pi < x < pi)的傅立葉級數F

10.範例11-5：求週期函數Square Wave方波函數f(x)，在(-pi < x < pi)的傅立葉級數Fourier Series


#重點1：要比較這題的標準解（以Σsigma的方式表示），與由 sypmy計算的差異
#重點2：如何在sympy裡面建立Square Wave方波函數
#重點3：了解傅立葉級數模擬函數會出現的一個現象：gibbs phenomenon 吉布斯跳動現象
#重點4：如何在sympy裡，消除傅立葉級數的gibbs phenomenon吉布斯跳動
# p = Piecewise((0, x < -1), (f, x <= 1), (g, True))
#本題的方波函數square wave
#f(x) = 1， 當 0 < x < pi
#f(x) = -1， 當 -pi < x < 0

from sympy import *
x = symbols('x')
f = Function('f')(x)
#重點2：如何在sympy裡面建立Square Wave方波函數：使用分段函數Piecewise()
#語法：Piecewise((x範圍1, x>0), (x範圍2, x<0))
#語法：局部區段的設定值：Piecewise((x範圍1, (0 f = Piecewise((1, (0 #plot(f)

#把方波函數轉成傅立葉級數（以(-pi < x < pi)週期）
s = fourier_series(f, (x, -pi, pi))
#只取n次級數
fs = s.truncate(n=15)
#重點3：了解傅立葉級數模擬函數會出現的一個現象：gibbs phenomenon 吉布斯跳動現象
#(1).傅立葉級數在f值不連續點(間斷處)，會出現所模擬的圖形大幅度波動現象，雖然把級數n提高，只能改善波動密度，但無法改善波動振幅
#(2).這種在不連續處的大幅度跳動，就是gibbs phenomenon 吉布斯跳動現象

#重點4：如何在sympy裡，消除傅立葉級數的gibbs phenomenon吉布斯跳動
#方法：把你的傅立葉級數，轉成sigma_approximation()，就能夠消除gibbs phenomenon吉布斯跳動
fssigma = s.sigma_approximation(15)

#畫出圖比較（fs vs fssigma）
p1 = plot(fs, fssigma, show=False)
p1[0].line_color = 'green'
p1[1].line_color = 'red'
p1[0].label = '$fourier(n=15)$'
p1[1].label = '$sigma_approximation(15)$'
p1.legend = True
p1.show()

#結果1：傅立葉級數的圖形，在斷點出，出現大幅度跳動(綠色圖)
#結果2：但是結果sigma_approximation轉換後，就可以消除gibbs phenomenon吉布斯跳動了



成果圖示：
程式碼內容

範例11-6：求週期函數f(x)=x**2，在(-pi < x < pi)的傅立葉級數

11.範例11-6：求週期函數f(x)=x**2，在(-pi < x < pi)的傅立葉級數Fourier Series


#注意：比較這題的標準解（以Σ:sigma的方式表示），與由 sypmy計算的差異
from sympy import *
x = symbols('x')
f = Function('f')(x)
f = x**2
#求出傅立葉級數
s = fourier_series(f, (x, -pi, pi))
#取出n=8的傅立葉級數
fs = s.truncate(n=8)
#轉成Σ:sigma _approximation
fssigma = s.sigma_approximation(8)
#印出數學表達式，與標準數學解來做比對
from sympy.interactive import printing
printing.init_printing(use_latex=True)
display(fssigma)

#印出圖形
p1 = plot(fs, fssigma, show=False)
p1[0].line_color = 'red'
p1[1].line_color = 'green'
p1[0].label = '$Fourier Series$'
p1[1].label = '$Σ:sigma_approximation$'
p1.legend = True
p1.show()


比較這題的標準解（以Σ:sigma的方式表示）

成果圖示：
程式碼內容

範例11-7：求週期2L函數f(x)=|x|，在(-3 < x <3)的傅立葉級數

12.範例11-7：求週期2L函數f(x)=|x|，在(-3 < x <3)的傅立葉級數Fourier Series

#物理很多函數是2L週期，不是2𝝿週期
from sympy import *
x = symbols('x')
f = Function('y')(x)
f = abs(x)
#傅立葉級數
s = fourier_series(f, (x, -3, 3))
fs = s.truncate(n=8)

#畫圖
plot(fs, line_color='red')
#顯示數學式
from sympy.interactive import printing
printing.init_printing(use_latex=True)
display(fs)

成果圖示：
程式碼內容

範例11-8：求週期Square Wave方波函數f(x)，在(-2 < x < 2)的傅立葉級數

13.範例11-8：求週期Square Wave方波函數f(x)，在(-2 < x < 2)的傅立葉級數Fourier Series


#本題的方波函數square wave
#f(x) = 0， 當 -2 < x < -1
#f(x) = 1， 當 -1< x < 1
#f(x) = 0， 當 1 < x < 2
from sympy import *
x = symbols('x')
f = Function('f')(x)
#注意：使用Piecewise,複合局部區域，1< x < 2，的表示方法為：(x<-1) & (x>-2)
f = Piecewise((0, (x>1)&(x<2)), (1, (x>-1)&(x<1)), (0, (x<-1)&(x>-2)))
#plot(f)
#轉成傅立葉級數
s = fourier_series(f, (x, -2, 2))
#截短成n=8, n=30
fs8 = s.truncate(n=8)
fs30 = s.truncate(n=30)
#畫圖
p1 = plot(fs8,fs30,show=False)
p1[0].line_color = 'red'
p1[1].line_color = 'green'
p1[0].label = '$fourier(n=8)$'
p1[1].label = '$fourier(n=30)$'
p1.legend=True
p1.show()

#結果1：不管n=8，或是n=30，都會有gibbs phenomenon吉布斯跳動

成果圖示：
程式碼內容

範例11-9：傅立葉積分：求週期為∞的方波函數f(x)，)的傅立葉積分

14.範例11-9：傅立葉積分：求週期為∞的方波函數f(x)，)的傅立葉積分Fourier Integration

#重點1：顯示的物理世界，很多波不是週期波，而是獨立波
#重點2：當不再說週期函數，而是單一波形函數時，傅立葉級數，就會被改成傅立葉積分（因為範圍∞，級數就會形成積分式）
#重點3：傅立葉積分的語法：The Mellin transform
#Mellin transform F(s) of f(x),
#F(s)=∫∞0xs−1f(x)dx. ...
#The Mellin transform is related via change of variables to the Fourier transfor

成果圖示：
成果圖示：
程式碼內容

1.何謂傅立葉變換

1.何謂傅立葉變換：
（1）對某非週期函數f(x)積分，但以形式表示

其中，自變數x表示時間（以秒為單位），
而變換後的變數ξ表示頻率（以赫茲為單位）

（2）傅立葉變換（Fourier transform）是一種線性積分變換，用於信號在時域（或空域）和頻域之間的變換
原函數坐標系為時域，水平軸為時間x
轉換後的坐標系為頻域，水平軸為頻率ξ（以赫茲為單位）



（3）轉換成頻譜函數
f：稱作原函數
f^： 經傅立葉變換生成的函數，稱為頻譜。

（4）可逆性
在許多情況下，傅立葉變換是可逆的，即可通過 f^，轉換得到其原函數f。

（5）轉換後的函數：是個複數函數（振幅 + 相位）
通常情況下，f 是實數函數，
而 轉換後的函數f^則是複數函數，用一個複數來表示振幅和相位

（6）物理意義：
A).也就是原本訊號波，是隨時間而變化的
經過傅立葉轉換後，可以看出這個波，是有很多不同頻率的波合成的，可以看出不同頻率波的權重值與振幅。

B).若是聲音的訊號，經過轉換後，若是出現高頻率的波，就可以判斷那是噪音，可以將之刪除。

（7）應用：
傅立葉變換在醫學、數據科學、物理學、聲學、光學、結構動力學、量子力學、數論、組合數學、機率論、統計學、訊號處理、密碼學、海洋學、通訊、金融等領域都有著廣泛的應用。
例如在訊號處理中，傅立葉變換的典型用途是將訊號分解成振幅分量和頻率分量。

2.傅立葉轉換公式有兩種

2.傅立葉轉換公式有兩種：

（1）傅立葉轉換公式(1)：以弧頻率ξ（赫茲）表示

其中，自變數x表示時間（以秒為單位），
而變換後的變數ξ表示頻率（以赫茲為單位）

（2）.傅立葉轉換公式(2)：以角頻率ω（弧度每秒）表示
傅立葉變換也可以寫成角頻率形式： ω = 2πξ其單位是弧度每秒


（3）原始函數f(t)的相位與頻率
如果原始波是：sin(θ)
不同頻率波：sin(θ/f)，例如：sin(θ)，sin(2θ)，sin(3θ)，sin(4θ)...就是不同頻率f的波
不同相位波：sin(θ+φ)，sin(θ+2φ)+sin(θ+3φ)....就是不同相位φ的波

3.傅立葉變換：就是兩個坐標的變換

3.傅立葉變換：就是兩個坐標的變換
（1）所謂的變換，就是兩個坐標的變換
（2）傅立葉認為：任一的週期函數，都可以用正弦餘弦函數（就是ex+iey)的指數複數函數）組合而成
（3）我們平時所看到的時域波形圖：f(t)
若是把它拉開，看它的所有組合波，
就會拉到另外一個坐標軸：g(ω），或 p(ξ)
其中：
弧頻率ξ（赫茲）
角頻率ω（弧度每秒）

示意圖：


（4）所以：傅立葉變換，就是坐標的轉換
由f(t) → A(ω）


（5）但轉換後的系統，有兩個坐標體系
g(ω) =（A(ω） + P(ω））= 頻譜 + 相位譜
A).一個坐標系統： A(ω）：可以看出不同頻率的『振幅』分量值關係（這個分量乃是此頻率的波與原波f(t)的相似度，越相似，值越大）
B).另外一個坐標系統：P(w)，不同頻率的『相位（phase）φ』
例如：
在角頻率ω=0處，可以算出C0，a0的值，這個就是初始相位φ
每個諧波的相位不一定相同（C0之不同）

4.Inverse Fourier Transform，反向傅立葉轉換

4.Inverse Fourier Transform，反向傅立葉轉換
（1）正向轉換：
由函數波的時域f(t)
→ 可以轉換成另外頻域坐標：g(ω) =（A(ω） + P(ω））= 頻譜 + 相位譜

（2）反向變轉換：
由頻域坐標：g(ω) =（A(ω） + P(ω））
→ 也可以反向變回時域f(t)的函數波

5.什麼時候要用到傅立葉變換？

5.什麼時候要用到傅立葉變換？
（1）當我們要分析其頻率的效果時：看頻譜分析
例如：某商店的整年份銷售數據，f(t)，t是每日，f是營收。
如果，我們想分析，每季，每月，每週，每日的營業額變化，該如何做？
方法：把數據轉成傅立葉變換→ FT = A(ω）
ω是把t換成頻率(季，月，週，日）

（2）需要做訊號處理時（例如：濾波，卷積，壓縮圖片，去噪）

6.傅立葉轉換的可視化

6.傅立葉轉換的可視化
（1）傅立葉轉換的可視化
https://www.youtube.com/watch?v=r18Gi8lSkfM
https://www.youtube.com/watch?v=spUNpyF58BY
https://www.youtube.com/watch?v=spUNpyF58BY&t=938s

7.傅立葉級數，傅立葉轉換的關聯與差異

7.傅立葉級數，傅立葉轉換的關聯與差異
（1）傅立葉級數：
有兩種方式可以表達傅立葉級數：正弦波餘弦波sin(x)+cos(x)，或是指數複數函數exp(iwt)
A).複雜的週期函數f，可以用一系列簡單的正弦波sin(x)、餘弦波(cos(x)的總和表示
自變數是x

B).以指數複數函數exp(iwt)來表示


（2）傅立葉轉換：


傅立葉變換將函數的時域（紅色）與頻域（藍色）相關聯。
頻譜中的不同成分頻率在頻域中以峰值形式表示

（3）為什麼需要傅立葉轉換：
因為：為了處理非週期性波，才需要傅立葉轉換

傅立葉分析最初是研究週期性現象，即傅立葉級數的，
後來通過傅立葉變換將其推廣到了非週期性現象。
理解這種推廣過程的一種方式是將非週期性現象視為週期性現象的一個特例，即其週期為無限長

(4)三種訊號坐標：
時域信號，
角頻率表示的傅立葉變換，
弧頻率表示的傅立葉變換

8.傅立葉變換家族

8.傅立葉變換家族：
下方列出了傅立葉變換家族的成員。容易發現，函數在時（頻）域的離散對應於其像函數在頻（時）域的週期性，反之連續則意味著在對應域的信號的非週期性

9.傅立葉級數(FS) vs 傅立葉轉換(FT)的應用差異

9.傅立葉級數(FS) vs 傅立葉轉換(FT)的應用差異：
（1）若訊號是週期訊號：則使用傅立葉級數（FS，Fourier Series)
（2）若訊號是非週期訊號：則使用傅立葉轉換（FT，Fourier Transform)



（3）若訊號是離散的非週期選號，則使用離散性訊號傅立葉轉換(D.T.F.T)
■注意：離散時間訊號常常可以用下圖表示

■注意：很多實際訊號是非週期，離散性的，所以D.T.F.T.常用

（3）轉換後的波，一個是離散的，一個是連續的
傅里葉級數，在時域是一個週期且連續的函數，而在頻域是一個非週期離散的函數。
傅里葉變換，將時域非週期的連續信號，轉換為一個在頻域非週期的連續信號。
(這是因為傅立葉變換，乃是對∞的積分，所以能夠包含所有的頻率，造成連續性的頻譜）

10.傅立葉積分 vs 傅立葉轉換的差異

10.傅立葉積分 vs 傅立葉轉換的差異
（1）相同點：
A).都是非週期訊號的處理。
B).或是說，是週期為∞
C).因為範圍為∞，所以用積分計算，不用級數計算
D).都是以

（2）差異點：
A).傅立葉積分：對非週期訊號方波函數f(x)的積分，但以正弦波餘弦波sin(x)+cos(x)來表示


B).傅立葉轉換：也是對非週期訊號方波函數f(x)的積分，但以指數複數函數exp(iwt)表示


（3）結論：這兩者的內容實際上述類似的。
週期為∞的函數 f(x) 之 Fourier 積分
也可以稱為：無限區間之 Fourier 轉換

（4）若非週期函數，只討論有限區間
若所討論的問題之定義域為有限區間，則可選用合適之有限區間的 Fourier 轉換加以解析；
若所討論的問題之定義域為無限區間，則可選用合適之無限區間的Fourier 轉換加以解析

11.波的表示式有兩種：歐拉公式：e^(ix) = cos(x) + i*sin(x)

11.波的表示式有兩種：歐拉公式：e^(ix) = cos(x) + i*sin(x)
（1）波的表示式有兩種：
A).以正弦波餘弦波sin(x)+cos(x)來表示
B).以指數複數函數exp(iwt)表示

（2）這兩者的關係式：歐拉公式
歐拉公式：e^(ix) = cos(x) + i*sin(x)
歐拉公式（Euler's formula，尤拉公式）是複分析領域的公式，它將三角函數與複指數函數關聯起來
e^{ix} = cos x + i*sin x

其中e 是自然對數的底數，i 是虛數單位，而 cos 和 sin 則是餘弦、正弦對應的三角函數，參數 x 則以弧度為單位
{x 為複數時仍然成立


（3）歐拉公式在數學、物理和工程領域應用廣泛

（4）當x為實數時，函數e^{ix}可在複數平面描述一單位圓。且x為此平面上一條連至原點的線與正實數軸的交角。

（5）歐拉公式在此提供複點至極坐標的變換:
任何複數z=x+yi皆可記為
z = x + iy = |z| (cos(phi) + i*sin(phi) = |z| e^{i*phi}
其中：x = Re{z}為實部。y = Im{z}為虛部



（6）e^(iωt)的順轉，與逆轉
e^(iωt)代表逆時針轉
e^(-iωt)代表順時針轉

（7）歐拉公式所描繪的，是一個隨著時間變化，在復平面上做圓周運動的點，隨著時間的改變，在時間軸上就成了一條螺旋線。
如果只看它的實數部分，也就是螺旋線在左側的投影，就是一個最基礎的餘弦函數。而右側的投影則是一個正弦函數

12.傅立葉轉換（FT）後是個複數： g(ω) = A(ω）+ iP(ω）

12.傅立葉轉換（FT）後是個複數： g(ω) = A(ω）+ iP(ω）= 頻譜 + i*相位譜
（1）一般所認為的FT = g(ω) = A(ω）
這個轉成頻域的坐標g(ω)，只是轉換後複數的實部（是個連續的波）

（2）FT，另外還有個虛部：i P(ω），表示相位
g(ω) = A(ω）+ iP(ω）= 頻譜 + i*相位譜

13.傅立葉轉換的應用：相機的美顏功能，與變聲

13.傅立葉轉換的應用：相機的美顏功能，與變聲
（1）把教室裡的聲音波，轉換成傅立葉變換→ g(ω)= A(ω）+ iP(ω）
其中的：
低頻波：男生說話
中高頻波：女生說話
高頻波：噪音

把高音頻刪除，再inverse逆轉回原波函數，就會把噪音消除。

（2）照片的美顏，或圖片的壓縮
一張2D空間的人之照片f(x=空間）→傅立葉轉換→ g(ω)= A(ω）+ iP(ω）
其中的：
低頻波：表示人的輪廓
高頻波：表示細節（斑點，青春痘）
若是在 g(ω)，把高頻訊號刪除，再轉回原函數f(x）→斑點就不見了

14.python有三種模組都支援快速傅立葉轉換

14.python有三種模組都支援快速傅立葉轉換（FFT,fast fourier transform)
（2）scipy：支援傅立葉級數，傅立葉轉換，快速傅立葉轉換FFT
優點：在處理離散訊號的FFT，支援很多函數，而且處理效率最快，建議使用

（2）numpy:支援傅立葉級數，傅立葉轉換，快速傅立葉轉換FFT
在處理離散訊號的FFT，支援很多函數

（3）sympy：支援傅立葉級數，傅立葉轉換，快速傅立葉轉換FFT
缺點1：在FFT支援的相關函數太少，很多離散性的訊號，無法處理，不建議使用
缺點2：訊號的點數只能處理少量數目（例如20個），若是超過太多（200個），就會當機無法計算

（4）scipy的傅立葉轉換範例教學：
https://docs.scipy.org/doc/scipy-0.17.0/reference/tutorial/fftpack.html
https://www.oreilly.com/library/view/elegant-scipy/9781491922927/ch04.html

範例12-1：求週期為∞的指數函數f(x)=exp(-x**2)的傅立葉轉換

21.範例12-1：求週期為∞的指數函數f(x)=exp(-x**2)的傅立葉轉換Fourier Transform
#觀念：現實的物理世界，很多波不是週期波，而是獨立波
#FT： 非週期性的連續波之傅立葉轉換語法： fourier_transform(f, x, k).doit()
#Inverse FT： 非週期性的連續波之傅立葉轉換語法： InverseFourierTransform(ft, k, x).doit()
from sympy import *
#from sympy.plotting import plot, PlotGrid
x = symbols('x')
k = symbols('k')
f = Function('f')(x)
f = exp(-x**2)

#(1).FT：傅立葉轉換
#fourier_transform(f, x, k)：表達式
#fourier_transform(f, x, k).doit()：把表達式simplify成簡約式，才能夠後續應用（例如畫圖，計算）
ft = fourier_transform(f, x, k).doit()
#print(ft)

#(2).inverse FT：反向傅立葉轉換成原始函數
#InverseFourierTransform(ft, k, x)：反向FT的表達式
#fourier_transform(f, x, k).doit()：把反向表達式simplify成簡約式，才能夠後續應用（例如畫圖，計算）
ftinverse = InverseFourierTransform(ft, k, x).doit()
#注意：若是沒有doit()，最後無法畫出圖，會出現錯誤

#(3).畫出三個子圖，比較：時域 → 頻域 → 反向時域
#sympy 畫出三個子圖的方法：使用plotting.PlotGrid(上下數目, 左右數目 , p1, p2, p3)
#sympy 顯示圖例legend的方法：若用label='FT=frequency spectrum',legend=True無法顯示
#必須另外寫p1[0].label = '$f(x)=exp(-x**2)$'，才能顯示出legend
#畫出時域：原始f(x)圖
p1 = plot(f, label='$f(x)=exp(-x**2)$', legend=True)
p1[0].label = '$f(x)=exp(-x**2)$'
#畫出頻域：FT頻譜圖
p2 = plot(ft,line_color='red', label='$FT=frequency spectrum$',legend=True)
p2[0].label = '$FT=frequency spectrum$'
#畫出反向時域
p3 = plot(ftinverse, line_color='green', label='$Inverse FT$',legend=True)
p3[0].label = '$Inverse FT$'
#在畫圖網格上分為上下三個圖
plotting.PlotGrid(3, 1 , p1, p2, p3)

成果圖示：
程式碼內容

範例12-2：D.T.F.T：離散資料，非週期數據，的傅立葉轉換法FFT

22.範例12-2：D.T.F.T：離散資料，非週期數據，的傅立葉轉換法FFT(Fast Fourier Transform)
#sympy要處理：離散性，非週期數據，乃是用快速傅立葉轉換法FFT(Fast Fourier Transform)
#語法：fft(seq, dps=None)
#注意：sympy要給定定的離散數據量，要是2*n倍，若不夠，系統會補0
#缺點：sympy的fft()函數，不好用。若給定太多的x數據，計算會變得很複雜，最後程式會當機
#使用sympy的FFT函數來做快速傅立葉轉換，有很多個嚴重缺點：
#（1）.讀入的離散資料，不能太多，讀入20個還可以計算，讀入200個數據就會無法計算
#（2）.fft轉換後的數據資料都是未簡化的數，若要將之轉成估值（numerical evaluation)會很麻煩
#（3）fft後續的陣列，sympy不提供對整個array陣列list的直接萃取數值（例如：只能對單一值做估值：.evalf()或是N(..)
#建議：使用：numpy或scipy，它們都有直接對陣列list操作的各種函數，而且轉換後的fft陣列，也是自動進行估值計算了

#Inverse FT： 非週期性的連續波之傅立葉轉換語法： InverseFourierTransform(ft, k, x).doit()
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
#from sympy.plotting import plot, PlotGrid
x = symbols('x')
k = symbols('k')
f = Function('f')(x)
#(1)原始量測數據
time_step = 0.1
time_vec = np.arange(0, 0.8, time_step)
print('(1).time=', time_vec)

sig = [0, 0.707, 1, 0.707, 0, -0.707, -1, -0.707]
#sig = [-0.2, 0.757, 0.95, 0.687, 0.02, -0.727, -0.95, -0.727]
print('(2).signal data =', sig)
plt.subplot(3,1,1)
plt.plot(time_vec, sig, label='raw data')
plt.legend()

#(2). FFTFFT：把時域，轉成頻域
#使用sympy的FFT快速傅立葉轉換函數fft，轉換後是個複數
sig_fft = fft(sig)
#print(sig_fft)
#這些複數都是未被簡化或估值的數據

#算出fft的水平頻率軸
#因為sympy這裡提供的函數不夠，所以使用scipy的函數
from scipy import fftpack
sample_freq = fftpack.fftfreq(len(sig), d = time_step)
print('sample_freq=',sample_freq)

#算出fft的垂直頻率軸振幅
#amplitude = np.abs(sig_fft)
#發現：np.abs()，有些複數沒有計算，所以必須要手動計算
#print('abs=', amplitude[0],sample_freq[0],len(amplitude),len(sample_freq))
#amplitude = np.sqrt(np.real(sig_fft)**2 + np.imag(sig_fft)**2)
amplitude = []
for item in sig_fft:
amplitude.append(N(Abs(N(item))))
#sympy數值估值函數：N()，或是.evalf()
#sympy計算複數絕對值 = sqrt(a**2+b**2) = Abs() = 振幅
print('amplitude=',amplitude)
plt.subplot(3,1,2)
plt.plot(sample_freq, amplitude, label='FFTFFT(frequency spectrum)')
plt.legend()

#(3).inver FFT：：把頻域，反轉成時域
sig_inverse_fft = ifft(sig_fft)
print('inverse fft=', sig_inverse_fft)
#這些複數都是未被簡化或估值的數據
plt.subplot(3,1,3)
plt.plot(time_vec, sig_inverse_fft,'r-', label='inverse fft data')

plt.legend()

成果圖示：
程式碼內容

範例12-3：D.T.F.T：使用spicy轉換離散資料，非週期數據的傅立葉轉換法FFT

23.spicy的#範例12-3：D.T.F.T：使用spicy轉換離散資料，非週期數據的傅立葉轉換法FFT(Fast Fourier Transform)
#reference:https://www.twblogs.net/a/5b8d116d2b717718833a60d8
#numpy與spicy都可以出來傅立葉轉換，而且提供很多支援的函數
#建議：因為spicy執行的效率較快，所以建議用spicy來處理快速傅立葉轉換
#(1).spicy的FFT語法：scipy.fftpack.fft(sig)
#或是：scipy.fftpack.fft(sig, n = 10)
#第一個參數sig：輸入的訊號序列或陣列list
#第二個參數n：輸入的訊號點數目 = len(sig) = sig.size
#注意：n最好是2*n倍，若不夠，系統會補0

#(2).把轉換後fft，由小到大排列，語法：scipy.fftpack.fftshift(sig_fft)
#注意：若用plot()畫frequency spectrum，則會自動大小排序，不需奧在fftshift()

#(3).計算頻譜圖的頻率陣列，語法：fftpack.fftfreq(n, d = time_step)
#其中 n = 訊號點數 = len(sig) = sig.size
#time step = d = Δt
#採樣點頻率 = fs = 1/d
#採樣點的頻率分辨率 = k*fs/N = k*1/(d*N)
#範例：sample_freq = fftpack.fftfreq(len(sig), d = time_step)

#(4).把頻率陣列重新排序(讓0頻率在中央)，語法：sample_freq_shift = fftpack.fftshift(sample_freq)
#不但頻率可以把0放在中央，fft後的複數函數也會依據頻率的排序而調整順序（會讓0頻率在中央）

#(5).計算頻譜圖的振幅陣列：因為頻率的順序已經改變，所以amplitude也要改版
#amplitude 公式：amplitude_shift = np.sqrt(np.real(sig_fft)**2 + np.imag(sig_fft)**2)

#(6).手動計算頻譜圖的振幅陣列(若是沒有shift頻率，那麼就可以自行手動計算amplitude)：
#方法二：Amplitude = np.abs(sig_fft)
#方法二：
# amplitude = []
# for item in sig_fft:
# amplitude.append(N(Abs(N(item))))
#sympy數值估值函數：N()，或是.evalf()
#sympy計算複數絕對值 = sqrt(a**2+b**2) = Abs() = 振幅

#(7).Inverse FT 語法： sig_inverse_fft = ifft(sig_fft)

from sympy import *
import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt
#x = symbols('x')

#(1)原始量測數據
time_step = 0.1
time_vec = np.arange(0, 0.8, time_step)
print('(1).time=', time_vec)

sig = [0, 0.707, 1, 0.707, 0, -0.707, -1, -0.707]
#sig = [-0.2, 0.757, 0.95, 0.687, 0.02, -0.727, -0.95, -0.727]
print('(2).signal data =', sig)
plt.subplot(3,1,1)
plt.plot(time_vec, sig, label='raw data')
plt.legend()

#(2). FFTFFT：把時域，轉成頻域
#使用sympy的FFT快速傅立葉轉換函數fft，轉換後是個複數
sig_fft = fftpack.fft(sig)
print('sig_fft=', sig_fft)
#依據頻率的排序而調整順序（會讓0頻率在中央）
sig_shift =fftpack.fftshift(sig_fft)
print('sig_shift=', sig_shift)

#算出fft的水平頻率軸
#因為sympy這裡提供的函數不夠，所以使用scipy的函數
sample_freq = fftpack.fftfreq(len(sig), d = time_step)
print('sample_freq=',sample_freq)
#把頻率陣列重新排序(讓0頻率在中央)，語法：fftpack.fftshift(sample_freq)
sample_freq_shift =fftpack.fftshift(sample_freq)
print('sample_freq_shift=',sample_freq_shift)

#算出fft的垂直頻率軸振幅
#amplitude = np.abs(sig_fft)
#發現：np.abs()，有些複數沒有計算，所以必須要手動計算
#print('abs=', amplitude[0],sample_freq[0],len(amplitude),len(sample_freq))

#另外，因為你把frequency的順序shift了，所以amplitude也必須依照sample_freq_shift的順序，手動計算
#amplitude 公式：amplitude_shift = np.sqrt(np.real(sig_fft)**2 + np.imag(sig_fft)**2)
amplitude_shift = np.sqrt(np.real(sig_shift)**2 + np.imag(sig_shift)**2)
#amplitude_shift = []
#for item in sig_shift:
# amplitude_shift.append(N(Abs(N(item))))
#sympy數值估值函數：N()，或是.evalf()
#sympy計算複數絕對值 = sqrt(a**2+b**2) = Abs() = 振幅
print('amplitude=',amplitude_shift)
plt.subplot(3,1,2)
plt.plot(sample_freq_shift, amplitude_shift, label='FFTFFT(frequency spectrum)')
plt.legend()

#(3).inver FFT：：把頻域，反轉成時域
#注意要invers，不可以針對sig_shift，因為它的順序以及被改變了
sig_inverse_fft = fftpack.ifft(sig_fft)
sig_inverse_fft_shift =fftpack.fftshift(sig_inverse_fft)
print('inverse fft=', sig_inverse_fft_shift)
#這些複數都是未被簡化或估值的數據
plt.subplot(3,1,3)
plt.plot(time_vec, sig_inverse_fft,'r-', label='inverse fft data')
plt.legend()

成果圖示：
成果圖示：
程式碼內容

範例12-4：把雜訊濾波：D.T.F.T：使用scipy做離散非週期性數據的傅立葉轉換

24.範例12-4：把雜訊濾波：D.T.F.T：使用scipy做離散非週期性數據的傅立葉轉換（FFT， Fast Fourier Transform)，有範例
#目的：我們給定一些以sin(x)為基礎的訊號，但隨機加上干擾噪音雜訊，當轉成頻譜後，把雜訊波消除，然後再inverse回原波
from sympy import *
import matplotlib.pyplot as plt
from scipy import fftpack
import numpy as np
#(1).每隔0.05秒測量一點，共200個數據
time_step = 0.05
time_vec = np.arange(0, 10, time_step)
#print(time_vec)

period = 5
#print(2*np.pi*time_vec/period)
#(2).以sin(2πt/p)為基礎的訊號，但隨機加上干擾噪音雜訊
sig = (np.sin(2*np.pi*time_vec/period) + 0.25*np.random.randn(time_vec.size))
#print(sig)
#print(np.round(sig,2))

#(3).把訊號用FFT轉成頻譜sig_fft複數
sig_fft = fftpack.fft(sig)
Amplitude = np.abs(sig_fft)
Power = Amplitude**2
Angle = np.angle(sig_fft)

#(4).計算出頻率軸(由sig_fft複數來計算)（不需要再shift重新排序了，因為plot時會自動調整）
sample_freq = fftpack.fftfreq(sig.size, d= time_step)
#print('sig.size=',sig.size,'len(sig)=',len(sig))
#print(Amplitude)
#print(sample_freq)

#(5).計算頻譜圖中的最大振幅處的頻率值
#原理：因為原訊號是sin(x)波，轉成頻域後，只會出現一個頻率的振幅，其它都是雜訊
#方法一：稍微複雜
Amp_Freq = np.array([Amplitude, sample_freq])
#Amp_Freq[0] = Amplitude
#Amp_Freq[1] = sample_freq
#numpy指令：argmax()：求出陣列最大值的所在索引index
#Amp_Freq[0,:].argmax()：表示在Amplitude裡面找出最大振幅的所在索引index
Amp_position = Amp_Freq[0,:].argmax()
#然後再到#Amp_Freq[1] （=sample_freq）裡面去查它對應的頻率值
peak_freq = Amp_Freq[1, Amp_position]

#方法二：更簡單更直觀的方法
#Amp_position = Amplitude.argmax()
#peak_freq = sample_freq[Amp_position]
print('Amp_position=', Amp_position)
print('peak_freq=', peak_freq)

#(6).濾波：把最大振幅波頻率以外的全部的波，其振幅都設定為0（沒有雜訊）
#最大漲幅波頻率=0.2， -0.2，所以f > 0.2 and f <-0.2的波，其sig_fft，都設定為0
high_freq_fft = sig_fft.copy()
#不需要迴圈，直接對陣列進行判別與設定：arrayb1[i > 5] =0
#目的：把sig_fft裡面，所有頻率
high_freq_fft[np.abs(sample_freq) > peak_freq] = 0
#print(high_freq_fft)

#(7).反向推算原函數（由已經濾波後的fft=high_freq_fft來反推）
filter_sig = fftpack.ifft(high_freq_fft)

#(8).畫圖: 點數據list
#用matplotlib可以直接畫處list數據
#在同一圖fig，畫出二條線的方法：最後再plt.show()即可
plt.subplot(2,1,1)
plt.plot(time_vec, sig, 'r-',label='raw data')
plt.legend()
#濾波後，畫圖
plt.subplot(2,1,1)
plt.plot(time_vec, filter_sig, 'b-',label='del noise')
plt.legend()
#plt.show();

#頻域畫圖
#在下方畫出sub plot的方法
plt.subplot(2,1,2)
plt.plot(sample_freq, Amplitude, 'r-',label='FFT(frequency spectrum')
plt.legend()
plt.show();

成果圖示：
成果圖示：
程式碼內容

範例12-5：從2個sin合成波訊號中濾除雜訊

25-範例12-5：從2個sin合成波訊號中濾除雜訊
#訊號的頻率=1 kHz，訊號間的總量測時間是=1.5 seconds.
#訊號的頻率(Sampling frequency)：Fs = 1000 Hz
#訊號的間隔時間= T = 1/Fs = 0.001
#訊號數目(Length of signal)= L = 1500
#時間矩陣(time_vector)= t = (0:L-1)*T
#模擬訊號：由兩個sin(x)波組合而成（50 Hz sinusoid of amplitude 0.7）+（120 Hz sinusoid of amplitude 1）
#合成波 S = 0.7*sin(2πt*50) + 1*sin(2πt*120)
#加上干擾雜訊： 2*random.ranint(時間矩陣點數的數目) = 2*random.ranint(time_vector.size)
#轉換成頻域，會發現不止50Hz, 120Hz，還有很多頻率的雜訊，要將之濾除

import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt
#(1).設定時間軸矩陣向量，頻率，訊號數目
Fs = 1000
T = 1/Fs
time_vec = np.arange(0, 1.5, T)
L = time_vec.size
print('L=',L, ', T=', T, ', Fs=', Fs)
print('time_vec=', time_vec)

#(2).建立2個sin的組合波（50 Hz sinusoid of amplitude 0.7）+（120 Hz sinusoid of amplitude 1）
#sig2 = 0.7*np.sin(2*np.pi*time_vec*50) + 1*np.sin(2*np.pi*time_vec*120)
sig2 = 0.7*np.sin(2*np.pi*time_vec*50) + 1*np.sin(2*np.pi*time_vec*120)
#np.random.randn(n)：建立一個一維n元素亂數矩陣, 內容是自然分佈，平均=0， 變化量=1(distribution of mean 0 and variance 1)
#np.random.randn(m, n)：建立一個一維mxn亂數矩陣
sig_nose = 2*np.random.randn(L)
sig = sig2 + sig_nose
#sig = sig2
print('sig=', sig)

#(3).畫出時域訊號圖(show the first 50 signal)
plt.subplot(3,1,1)
plt.plot(1000*time_vec[0:50], sig[0:50], label='raw signal')
plt.legend()
#(4).用FFT=轉成頻域(用scipy的fft函數，效率較好)
sig_fft = fftpack.fft(sig)

#(5).計算頻域的水平軸頻率矩陣
#scipy語法：fftpack.fftfreq(sig_fft.size, d=訊號間隔時間)
sample_freq = fftpack.fftfreq(sig_fft.size, d=T)

#(6).計算頻域的振幅矩陣
#amp2 = np.array()
amplitude = np.abs(sig_fft)/L
#只看頻率為正值部分，振幅正負部位合在一起，2倍
amplitude2 = 2*amplitude

#(7).畫圖：頻域頻譜圖
plt.subplot(3,1,2)
#頻率報復只顯示一半:0～L/2，這是正值部分
plt.plot(sample_freq[0:int(L/2)], amplitude2[0:int(L/2)], color='red',label='fft(frequency)')
plt.legend()

#(8).印出振幅最大1st的頻率
index_peak_1stmax = amplitude2[:].argmax()
#print('index_peak_1stmax=', index_peak_1stmax)
freq_peak_1smax = sample_freq[index_peak_1stmax]
print('freq_peak_1smax=', freq_peak_1smax)
amplitude_peak_1smax = amplitude2[index_peak_1stmax]
print('amplitude_peak_1smax=', amplitude_peak_1smax)

#(9).濾波：把最大振幅波頻率以外的全部的波，其振幅都設定為0（沒有雜訊）
#最大漲幅波頻率, A < 0.5的波，其sig_fft，都設定為0
high_sig_ft = sig_fft.copy()
high_sig_ft[np.abs(amplitude2) < 0.4] = 0
#print(high_freq_fft)

#(10).反向Inverse推算原函數（由已經濾波後的fft=high_freq_fft來反推）
sig_inverse_fft = fftpack.ifft(high_sig_ft)

#(11).在下方畫出sub plot的方法
plt.subplot(3,1,3)
plt.plot(time_vec[0:50], sig_inverse_fft[0:50],'r-', label='inverse fft data')
plt.legend()
plt.show()

成果圖示：
成果圖示：
程式碼內容

檔案

1.日期時間函式庫

1.日期時間函式庫
(1).時間函式庫有三個：datetime，date，time模組
模組datetime：可以處理日期與時間函式庫，可以傳回目前的日期，時間
模組date：可以傳回目前的日期

範例13-1:模組datetime（顯示日期，時間

(2).範例13-1:模組datetime（顯示日期，時間）
公式：
from datetime import datetime
a1 = datetime.today()
y1 = a1.year
範例13-1.py
from datetime import datetime
a1 = datetime.today()
print('現在日期時間 =',a1)
y1 = a1.year
m1 = a1.month
d1 = a1.day
h1 = a1.hour
min1 = a1.minute
s1 = a1.second
print(y1,'年',m1,'月',d1,'日',h1,'時',min1,'分',s1,'秒')
程式碼內容

範例13-2:模組date（只能顯示日期）

(3).模組date（只能顯示日期）
公式：
from datetime import date
a1 = datetime.today()
y1 = a1.year

範例13-2.py
from datetime import date
a1 = date.today()
print('現在日期=',a1)
y1 = a1.year
m1 = a1.month
d1 = a1.day
print(y1,'年',m1,'月',d1,'日')
程式碼內容

範例13-3:模組timedelta（間隔日期，時間）

(4).模組timedelta（間隔日期，時間）
公式：
from datetime import datetime,timedelta
a1 = datetime.today()
delta2 = timedelta(days=5)
a2 = a1 + delta2

範例13-3.py
from datetime import datetime,timedelta
a1 = datetime.today()
print('今天是=',a1)
delta2 = timedelta(days=5)
a2 = a1 + delta2
print('5天後=',a2)
delta3 = timedelta(hours=3)
a3 = a1 +delta3
print('3小時後=',a3)
程式碼內容

範例13-4:傳回自1970年1月1日凌晨0：0：0開始至今的秒數

(5).模組time（時間）
公式：
import time
t1 = time.time()

功能：傳回自1970年1月1日凌晨0：0：0開始至今的秒數

範例13-4.py
import time
t1 = time.time()
print('1970年至今的秒數=',t1)
程式碼內容

範例13-5：計算呼叫某個函式所需要的時間

範例13-5.py
目的：計算呼叫某個函式所需要的時間
import time
def func():
[i for i in range(9999999)]
t1 = time.time()
func()
t2 = time.time()
print('函式執行時間 = ', t2-t1)
程式碼內容

2.可以迭代（iterate）的函式庫= itertools，函式enumerate，zip，filter，map，reduce

2.可以迭代（iterate）的函式庫
(1).種類：迭代模組或函式 = itertools，函式enumerate，zip，filter，map，reduce

(2).enumerate(a)是列舉函數（編號i，元素x）
☎功能：可以把每一個元素加上編號
☎配合list：list(enumerate(a)))
☎顯示一維串列：print((x for i,x in enumerate(a)))
☎顯示二維串列：print((x[0]+'姓名是'+x[1] for i,x in enumerate(a)))

範例13-6:顯示陣列

☎範例13-6.py
a=['tom','mike','peter','yellow']
#for顯示陣列
for i in a:
print(i)
#enumrate顯示陣列 enumerate
print([x for i,x in enumerate(a)])

#enumrate顯示陣列 enumerate
print(list(enumerate(a)))

b = (('tom','湯姆'),('mike','麥克'),('peter','彼德'))
#for顯示陣列
for i in b:
print(i)
for i in b:
print(i[0],',',i[1])
#enumrate顯示二維陣列 enumerate
print([x[0] for i,x in enumerate(b)])
程式碼內容

範例13-7:zip(a,b)是結合兩個可迭代物

(3).zip(a,b)是結合兩個可迭代物件
☎功能：結合兩個可迭代物件，對應結合成一個新的物件
☎必須配合list：list(zip(a, b)))

☎範例13-7.py
a = ['tom','mike','peter','yellow']
b = [95,85,66,75]
#結合兩個可迭代物件
c = zip(a, b)
print(list(c))
程式碼內容

範例13-8:模組itertools:給定起始值，遞增值，就可以產生無窮數列

(4).模組itertools
功能：給定起始值，遞增值，就可以產生無窮數列
☎範例13-8.py
import itertools
# itertools.count(起始值,遞增值)
nums = itertools.count(1,2)
#print(list(nums))
#結果： 1 3 5 7 9....

# itertools.cycle(資料列)
nums = [1,2,3]
newnums = itertools.cycle(nums)
#print(list(newnums))
#結果： 1 2 3 1 2 3 1 2 3....

# itertools.repeat(資料列,次數)
nums = [1,2,3]
newnums = itertools.repeat(nums,2)
#print(list(newnums))
#結果： [[1,2,3],[1,2,3],]....

# itertools.accumulate(資料列)
nums = [1,2,3,4,5,6]
newnums = itertools.accumulate(nums)
print(list(newnums))
#結果： [1, 3, 6, 10, 15, 21]

# itertools.chain(資料列)
newnums = itertools.chain('how','are you')
print(list(newnums))
#結果： ['h', 'o', 'w', 'a', 'r', 'e', ' ', 'y', 'o', 'u']

# itertools.combinations(資料列',取n個)
newnums = itertools.combinations('ABC',2)
print(list(newnums))
#結果： [('A', 'B'), ('A', 'C'), ('B', 'C')]

# itertools.permutations(資料列',取n個)
newnums = itertools.permutations('ABC',2)
print(list(newnums))
#結果： [('A', 'B'), ('A', 'C'), ('B', 'C')]
程式碼內容

範例13-9:filter(過濾函數def或lambda,資料列)

(5).過濾模組filter
公式：filter(過濾函數def或lambda,資料列)
功能：
☎範例13-9.py
#顯示偶數
a = [1,2,3,4,5,6,7,8,9,10]
b = filter(lambda x:x%2==0, a)
print(list(b))

#顯示有'o'的元素
a = ['tom','girl','man','took']
b = filter(lambda x:x.find('o')>=0, a)
print(list(b))
程式碼內容

範例13-10:過濾模組map

(6).過濾模組map
公式：map(過濾函數def或lambda,資料列)
功能：
☎範例13-10.py
#顯示偶數
def getnum(x):
if x%2 ==0:
return x a = [1,2,3,4,5,6,7,8,9,10]
b = filter(getnum, a)
print(list(b))

#顯示有'o'的元素
def getstr(x):
if x.find('o')>=0:
return x
a = ['tom','girl','man','took']
b = filter(getstr, a)
print(list(b))
程式碼內容

範例13-11:操作每個元素模組reduce

(7).操作每個元素模組reduce
公式：reduce(過濾函數def或lambda,資料列)
功能：逐一操作每個元素
☎範例13-11.py
from functools import reduce
#計算總和
def getsum(x,y):
return x+y
a = [1,2,3,4,5]
b = reduce(getsum, a)
print(b)

#連接所有的元素
a = ['t','r','i','g','e','r']
b = reduce(lambda x,y:x+y, a)
print(b)
程式碼內容

範例13-12:套件collections-模組OderedDict：可將兩個串列，合併形成成一個字典dict

3.套件collections
import collections
(1).模組OderedDict 功能：將兩個串列，合併形成成一個字典dict
☎比較：合併串列建立字典的兩種方法
☎範例13-12.py
import collections
a1 = ['姓名','住址','性別','電話']
b1 = ['張三','台南市忠孝東路15好','男','0912751477']
#合併串列建立字典，方法一
d1 = collections.OrderedDict(zip(a1,b1))
print(d1)
print('姓名=',d1['姓名'])

#合併串列建立字典，方法二：zip()合併兩個串列，轉成dict()
c1 = zip(a1,b1)
d1 = dict(c1)
print(d1)
print('姓名=',d1['姓名'],'電話=',d1['電話'])
程式碼內容

範例13-13:模組deque：操作可迭代資料列，新增，刪除，反轉

(2).模組deque
功能：操作可迭代資料列，新增，刪除，反轉
☎範例13-13.py
from collections import deque
a1 = [1,2,3,4,5]
#建立deque物件
d1 = deque(a1)
#反轉資料列
d1.reverse()
print('反轉=',d1)
#新增到資料列右邊 append(元素)
d1.append('R')
print('右邊新增=',d1)
#新增到資料列左邊 appendleft(元素)
d1.appendleft('L')
print('左邊新增=',d1)
#刪除資料列右邊一個 pop()
d1.pop()
print('刪除右邊1個=',d1)
#刪除資料列左邊一個 pop()
d1.popleft()
print('刪除左邊1個=',d1)
#刪除某個元素 remove(元素)
d1.remove(3)
print('刪除某個元素3=',d1)
d1.append(2)
#計算某個元素2出現的次數 count(元素)
n1 = d1.count(2)
print('計算某個元素2出現的次數=',n1)
程式碼內容

例13-14:模組dCounter：可計算可迭代資料列相同元素的次

(3).模組dCounter
功能：計算可迭代資料列相同元素的次數
☎範例13-14.py
from collections import Counter
s = 'today is a good day'
a1 = list(s)
print(a1)
#建立counter物件（字典）
c1 = Counter(a1)
print(c1)
#顯示所有『鍵』
print(c1.keys())
#顯示所有『值』
print(c1.values())
#顯示所有『元素』
print(list(c1.elements()))
#顯示所有『鍵值（tuple）』
print(c1.most_common())
程式碼內容

1.存取xml:模組：xml.etree.ElementTree

1.python的數學函數庫：math，sympy，cmath （1）math： 功能：乘方、開方、對數，冪函數與對數函數，三角函數，角度轉換，雙曲函數，特殊函數，常量，cos,sin,e,log,tan,pow Math模組 sqrt(x) 傳回 x 的平方根 pow(x, y) 傳回 x 的 y 次方 exp(x) 傳回 x 的自然指數 expm1(x) 傳回 x 的自然指數-1 (在 x 接近 0 時仍有精確值) log(x [, b]) 傳回 x 以 b 為基底的對數 (預設 b=e 自然對數) log10(x) 傳回 x 的常用對數 (以 10 為底數) degrees(x) 傳回弧度 x 的角度 (degree) radians(x) 傳回角度 x 的弧度 (radian) dist(p, q) 傳回兩個座標點 p, q 的歐幾里得距離 (畢式定理斜邊) hypot(coor) 傳回座標序列 coor 的歐幾里得距離 sin(x) 傳回 x 的正弦值 cos(x) 傳回 x 的餘弦值 tan(x) 傳回 x 的正切值 asin(x) 傳回 x 的反正弦值 (sin 的反函數) acos(x) 傳回 x 的反餘弦值 (cos 的反函數) atan(x) 傳回 x 的反正切值 (tan 的反函數) atan2(y, x) 傳回 y/x 的反正切值 (tan 的反函數)=atan(y/x) sinh(x) 傳回 x 的雙曲正弦值 cosh(x) 傳回 x 的雙曲餘弦值 tanh(x) 傳回 x 的雙曲正切值 asinh(x) 傳回 x 的反雙曲正弦值=log(x+sqrt(x**2+1)) acosh(x) 傳回 x 的反雙曲餘弦值=log(x+sqrt(x**2-1)) atanh(x) 傳回 x 的反雙曲正切值=1/2*log((1+x)/(1-x)) fabs(x) 傳回 x 的絕對值 (或稱模數, modulus) floor(x) 傳回浮點數 x 的向下取整數 (即小於 x 之最大整數) ceil(x) 傳回浮點數 x 的向上取整數 (即大於 x 之最小整數) trunc(x) 傳回浮點數 x 的整數部分 (捨去小數) modf(x) 傳回浮點數 x 的 (小數, 整數) 元組 factorial(x) 傳回 x 階乘 (x!, x=整數) gcd(x, y) 傳回整數 x, y 之最大公因數 comb(n, k) 傳回 n 取 k 的組合數 (不依序不重複) perm(n, k) 傳回 n 取 k 的組合數 (依序不重複) modf(x, y) 傳回 x/y 之精確餘數 (浮點數 float) fsum(iter) 傳回可迭代數值 iter 之精確總和 isclose(x, y) 若 a, b 值很接近傳回 True (預設差小於 1e-9) isfinite(x) 若 x 不是 nan 或 inf 傳回 True, 否則 False isnan(x) 若 x 為 nan 傳回 True, 否則 False isinf(x) 若 x 為 inf 傳回 True, 否則 False （2）sympy： Sympy是一個數學符號庫（sym代表了symbol，符號），包括了積分，微分方程，三角等各種數學運算方法，是工科最基本的數學函數庫，用起來媲美matlab，而且其精度比math函數庫精確。 功能： simplify運算式化簡，solve方程自動求解 limit求極限，diff求導 dsolve()計算微分方程 intergrate積分計算：1.定積分，2.不定積分，3.雙重定積分，4. 雙重不定積分 （3）cmath：專門用來處理複數運算。

範例14-1:讀取person.xml的所有每個節點資訊，查詢所有的mail，查詢卓水信資料

☎方法1:範例14-1.py
下載person.xml
#目的讀取person.xml的所有每個節點資訊，查詢所有的mail，查詢卓水信的個人所有資料
程式碼：
import xml.etree.ElementTree as xml
tree = xml.ElementTree(file='person.xml')
root = tree.getroot()
print('根目錄的標籤名稱',root.tag)

# root.tag 讀取 tag 名
# root.attrib 讀取 attributes
# root[0][1].rank 讀取文本值
#輸出根節點元素的tag名：person
print(root.tag)
#輸出根節點元素的attributes(空值)：{}
print(root.attrib)
#輸出第一個的三個子元素的值：
print(root[0][0].text)
print(root[0][1].text)
print(root[0][2].text)

#查找所有子孫元素：tag標籤 = name
for elem in root.iter('name'):
print(elem.tag, elem.attrib, elem.text)

#比較find() 和 findall()
#查找第一個子元素：find()
print(list(root.find('student')))
#查找所有子元素：findall()
print(list(root.findall('student')))

#查找所有子元素：findall()
#tree.findall('student/*') //查詢孫子節點元素
for elem in root.findall('student/*'):
print(elem.tag, elem.text)

#tree.findall('.//name') //查詢任意層次元素
for elem in root.findall('.//tel'):
print(elem.tag, elem.text)

#查詢 包含name屬性的student
#查詢tree.findall('student[@name]')
for elem in root.findall('student[@name]'):
print(elem.tag, elem[1].text)

#查詢 name屬性為卓水信的student
#tree.findall('student[name="卓水信"]')
for elem in root.findall('student[name="卓水信"]'):
print('id=',elem[0].text)
print('name=',elem[1].text)
print('tel=',elem[2].text)

#顯示 第一個student
for elem in tree.findall('student[1]'):
print(elem[0].text,elem[1].text,elem[2].text)

#顯示 第二個student
for elem in tree.findall('student[2]'):
print(elem[0].text,elem[1].text,elem[2].text)

#顯示最後一個student
#tree.findall('student[last()]')
for elem in tree.findall('student[last()]'):
print(elem[0].text,elem[1].text,elem[2].text)

#顯示倒數第二個country
#tree.findall('country[last()-1]')
for elem in tree.findall('student[last()-1]'):
print(elem[0].text,elem[1].text,elem[2].text)

#顯示第1階層所有階層的標籤名稱，屬性
print('顯示第1階層所有階層的標籤名稱，屬性：')
for child in root:
print(child.tag, child.attrib, child.text)

#顯示第2階層的標籤名稱，文字值
print('顯示第2階層的標籤名稱，文字值：')
for i,child1 in enumerate(root):
print(child1.tag,i)
for child2 in child1:
print(child2.tag, '=', child2.text)

#查詢某個tag標籤名稱
print('查詢某個tag標籤名稱：')
for elem in root.iterfind('student/mail'):
print(elem.tag,'=',elem.text)

#查詢某個tag標籤tel的所有文字
print('查詢某個tag標籤的所有文字：')
for elem in root:
print(elem.find('tel').text)

#查詢卓水信是否存在
print('查詢卓水信是否存在：方法1')
#方法1
for elem in root.findall('student[name="卓水信"]'):
if elem[1].text=='卓水信':
print(elem[1].text,'存在')
print('查詢卓水信是否存在：方法2')
#方法2
for elem in tree.iter():
if elem.text=='卓水信':
print(elem.tag,'=',elem.text,'存在')

#查詢卓水信的個人所有資料
#方法1
print('\n查詢卓水信的個人所有資料：方法1')
for elem in root.findall('student[name="卓水信"]'):
print(elem[0].text, elem[1].text,elem[2].text)
#方法2
print('\n查詢卓水信的個人所有資料：方法2')
for elem1 in tree.iter():
for elem2 in elem1:
if elem2.text=='卓水信':
print('name =',elem1.find('name').text)
print('tel =',elem1.find('tel').text)
print('mail =',elem1.find('mail').text)
程式碼內容

範例14-2:讀取person.xml的所有每個節點資訊，查詢所有的mail，查詢卓水信

☎方法2範例14-2.py
#目的讀取person.xml的所有每個節點資訊，查詢所有的mail，查詢卓水信的個人所有資料
程式碼：
import xml.etree.ElementTree as xml
tree = xml.ElementTree(file='person.xml')
#顯示所有的tag name = tree.iter()
print('顯示所有的tag name')
for elem in tree.iter():
print(elem.tag, elem.attrib,elem.text)
#尋找某個tag name的資料 = tree.iterfind('tag name')
print('尋找某個tag name')
for elem in tree.iterfind('student/mail'):
print(elem.tag,'=',elem.text)
#尋找某個tag text的資料 = tree.iterfind('tag name')
print('尋找某個tag attribute')
for elem in tree.iterfind('student[@hash="1cdf045c1"]'):
print(elem.tag,'=',elem.attrib)

print('\n尋找全部的mail')
for elem in tree.iterfind('student/mail'):
print(elem.tag,'=',elem.text)

#查詢卓水信的個人所有資料
print('\n查詢卓水信的個人所有資料：')
for elem1 in tree.iter():
for elem2 in elem1:
if elem2.text=='卓水信':
print('name =',elem1.find('name').text)
print('tel =',elem1.find('tel').text)
print('mail =',elem1.find('mail').text)
程式碼內容

範例14-3:修改並存入xml

(3).修改並存入xml文檔
☎範例14-3.py
#修改卓水信的個人id_no = 999999
程式碼：
import xml.etree.ElementTree as xml
tree = xml.ElementTree(file='person.xml')
root = tree.getroot()
#修改
print('\n修改卓水信的個人資料：')
for elem in root.findall('student[name="卓水信"]'):
elem[0].text = '999999'

#存檔xml
tree.write('person2.xml') #儲存
#顯示
for elem1 in root:
for elem2 in elem1:
print(elem2.tag,elem2.attrib, elem2.text)
程式碼內容

線上XML/JSON互相轉換工具

線上XML/JSON互相轉換工具：
http://tools.itread01.com/code/xmljson

線上格式化XML/線上壓縮XML：
http://tools.itread01.com/code/xmlformat

XML線上壓縮/格式化工具：
http://tools.itread01.com/code/xml_format_compress

XML程式碼線上格式化美化工具：
http://tools.itread01.com/code/xmlcodeformat

更多關於Python相關內容感興趣的讀者可檢視：
《Python操作xml資料技巧總結》、《Python資料結構與演算法教程》、《Python Socket程式設計技巧總結》、《Python函式使用技巧總結》、《Python字串操作技巧彙總》、《Python入門與進階經典教程》及《Python檔案與目錄操作技巧彙總》
https://www.itread01.com/article/1535600015.html

2.讀取網頁：request(url)

2.讀取網頁：request(url)
☎模組：import urllib.request as request
☎讀取網頁：web = request.urlopen(網址)
☎讀取網頁內容並解碼 = web.read().decode()

範例14-4:讀取網頁：web = request.urlopen(網址)

☎範例14-4.py
程式碼：
import json
import urllib.request as request
url = 'http://web.tsu.edu.tw/bin/home.php'
web = request.urlopen(url)
print('網址=',web.geturl())
print('狀態（200表示OK)=',web.status)
print('取得網頁標頭=',web.getheaders())
print('取得網頁標頭=',web.getheaders())
txt = web.read()
print('取得網頁內容(byte格式)=',txt)
print('取得網頁內容(解碼byte格式)=',txt.decode())
程式碼內容

3.存取 json(模組：json)

3.存取 json
(1).模組：json
import json

(2).把字典dict轉換成json
公式：json = json.dumps(字典)

(3).把json轉換成字典
公式：字典 = json.loads(json)

範例14-5:轉成json：jumps。轉成dict：loads

(4).☎範例14-5.py
程式碼：
import json
#把字典dict轉成json：jumps
a1 = {'tom':'0912456789','mike':'0965258741','peter':'0965789365'}
j1 = json.dumps(a1)
print(j1)

#把串列list轉成json：jumps
a2 = [['tom','湯姆'],['mike','麥克'],['peter','彼德']]
j2 = json.dumps(dict(a2))
print(j1)

#把json轉成字典dict：loads
#轉成上面的json = j1
a1 = json.loads(j1)
print('tom的電話 = ',a1['tom'])

#把json轉成字典dict：loads
#轉成上面的json = j2
a2 = json.loads(j2)
print(a2)
print('mike的電話 = ',a2['mike'])
程式碼內容

範例14-6：讀取網絡上的json檔案

(5).讀取網絡上的json檔案
☎範例14-6.py #目的：讀取網路上一個json檔案：http://acupun.site/lecture/jquery_phoneGap/json/book.json
程式碼：
import json
import urllib.request as request
url = 'http://acupun.site/lecture/jquery_phoneGap/json/book.json'
#開啟網頁：request.urlopen(url)
web = request.urlopen(url)
#讀取網頁文字：web.read()
txt = web.read()
# txt = web.read().decode()
#json轉成字典
dic1 = json.loads(txt)
print(dic1)
print('第一本書 = ',dic1[0]['title'], dic1[0]['author'])
程式碼內容

範例14-7：讀取電腦上的json檔案

(6).讀取電腦上的json檔案
☎範例14-7.py
程式碼：
import json
import os
#讀檔案
f1 = open('school.json','rt',encoding='utf-8-sig')
#讀檔案內容
txt = f1.read()
#print(txt)
#把json轉成字典dict
dict1 = json.loads(txt)
#print(dict1)
print('第一間學校 = ',dict1[0]['name'], dict1[0]['address'])
for elem in dict1:
print(elem['name'],elem['address'])
程式碼內容

3.網路爬蟲BeautifulSoup

3.網路爬蟲BeautifulSoup: 讀取並分析html網頁標籤
(1).先安裝第三方函數庫，使用：pip install beautifulsoup4
下載並安裝套件

(2).讀取下載在電腦的網頁
☎#注意：這個網頁.htm，必須放在電腦檔案內，不可讀取網絡上網址的網頁
☎（正確）fin = open('web1.htm',encoding='utf-8')
☎（錯誤）fin = open('https://www.python.org/',encoding='utf-8')
☎公式
from bs4 import BeautifulSoup as soup
fin = open('網頁.htm',encoding='utf-8')
txt = fin.read()
htm = soup（txt,'html.parse')

☎讀取網頁標題
print(htm.title.prettify())

☎讀取網頁標籤
公式：BeautifulSoup.find_all(tag)
找出標籤為tag的所有元素
例如：for item in htm.find_all('tr'):
print(item)

☎讀取網頁標籤的innerHtml文字：三種方法：
for item in htm.find_all('a',href='http://epaper.edu.tw/'):
print('innerhtml內容=',item.contents)
print('innerhtml內容=',item.contents[0])
print('innerhtml內容=',item.string)

☎讀取網頁標籤
公式：BeautifulSoup.find_all(tag, attr)
找出標籤為tag＋屬性為attr的所有元素
例如：for item in htm.find_all('tr',class='table_head'):
print(item)

範例14-8：讀取網頁標題

☎範例14-8.py
☎（正確）fin = open('web1.htm',encoding='utf-8')
☎（錯誤）fin = open('https://www.python.org/',encoding='utf-8')
程式碼：
from bs4 import BeautifulSoup as soup
fin = open('web1.htm',encoding='utf-8')
#fin = open('https://www.python.org/',encoding='utf-8')
txt = fin.read()
htm = soup(txt,'html.parser')
#讀取網頁標題
print(htm.title.prettify())
程式碼內容

範例14-9：讀取網址的網頁

(4).讀取網址的網頁：
☎範例14-9.py
#目的：讀取www.tsu.edu.tw網頁，如何顯示所有超連結的元素
#方法：找出標籤為tag的所有元素
例如：for item in htm.find_all('a'):
print(item)
程式碼：
import urllib.request as request
from bs4 import BeautifulSoup as soup
#使用request讀取網頁內容
url = 'http://web.tsu.edu.tw/bin/home.php'
web = request.urlopen(url)
txt = web.read().decode()
#print(txt)
#使用beautifulSoup發現網頁內容
htm = soup(txt, 'html.parser')
#讀取網頁標題
print(htm.title.prettify())
#讀取網頁所有的超連結
for item in htm.find_all('a'):
print(item)
#讀取網頁所有的超連結中，屬性是href='http://epaper.edu.tw/'
for item in htm.find_all('a',href='http://epaper.edu.tw/'):
print('指定屬性超連結a = ',item)
print('innerhtml內容=',item.contents)
print('innerhtml內容=',item.contents[0])
print('innerhtml內容=',item.string)
程式碼內容