陳擎文教學網:python求解數學式(高中數學,大學數學,工程數學,微積分)
方程式,多項式,函數,三角函數,微分,積分,泰勒級數,傅立葉級數,傅立葉積分,傅立葉轉換,一階微分方程式,二階微分方程式,三階微分方程式,偏微分方程式,聯立方程式矩陣解法,張量,向量的坐標轉換, 複數,極坐標

資源(Resource)

chp1.前言,執行python三種方法,安裝anaconda

chp2.python的數學函數庫:math,sympy,cmath,scipy,numpy

chp3.多項式,方程式,函數

chp4.複數complex number,與極坐標系統(Polar coordinate system)

chp5.微分,函數的微分

chp6.積分,integral

chp7.微分方程式,ODE:ordinary differential equation

chp8.二階常係數微分方程式,Second-Order Differential Equations

chp9.矩陣算法,解聯立方程式,線性代數

chp10.泰勒級數,泰勒多項式,用泰勒展開式模擬函數,Taylor Series

chp11.傅立葉級數,Fourier Series

chp12.傅立葉轉換(Fourier Transform, FT)

chp13.小波分析(wavelet analysis),小波轉換(wavelet transform)

chp14.偏微分方程式,PDE:partial differential equation


python資源
官網 python官網 vscode官網    
python 教學網站 w3c school(中文版) w3c school(英文版)    
線上執行python online

https://www.python.org/shell/(建議用這個)

https://repl.it/languages/python3


https://www.onlinegdb.com/online_python_compiler

https://www.tutorialspoint.com/execute_python_online.php

其它資源
  上課程式碼即時貼網站 線上黑板( Online blackboard) 上課即時貼(guestBook) 4月30日範例 有道翻譯
  考試題目 考試題目(Exam) 畢業門檻(Graduation threshold) 2018通識課規定 DWCS3
  Goole 輸入法(Input software) Goole輸入法(Input:exe) Goole 輸入法(Input:zip) online goole input(中文) online goole input(英文)
  Windows+Apache
+MySQL+PHP整合安裝系統
phpstudy2014 phpstudy2018 XAMPP wampserver
  免費的網頁編輯軟體 微軟的vs code(visual studio code)      
 

 
chp1.前言,安裝anaconda
1.使用python的三種方法 2.Anaconda下載點 3.安裝anaconda 4.Anaconda cmd指令
5.建立Anaconda虛擬環境 6.使用Spyter編譯器 7.網頁版python編輯器jupyter notebook 8.其它線上雲端可編譯的python平台

1.前言

Python堪稱是大數據與AI時代的最重要程式語言,在資料處理上有著非常重要的地位。而隨著AI的興起,讓傳統的零售業、金融業、製造業、旅遊業,以及政府都爭相投入,無不希望能運用數據分析與預測來協助決策方向,也讓新興的數據分析師、資料分析師成為熱門職業,因此本課程將講解如何使用網絡爬蟲技術以掌握資料爬取分析、視覺化呈現,以及儲存交換應用的關鍵技術。

Python資料處理的三大技術分別是:擷取分析、視覺化呈現與儲存應用。
而其應用的範疇包括:網路爬蟲、資料正規化、資料視覺化、資料儲存與讀取(CSV、Excel、Google試算表、SQLite、MySQL)、批次檔案下載、公開資料應用、API建立、驗證碼辨識。

Python大數據分析最重要的四個模組

1.Python大數據分析最重要的四個模組
Python資料分析最重要的四個模組:pandas、numpy、scipy、matplotlib。

(1)pandas:是基於numpy的資料分析工具,能夠快速的處理結構化資料的大量資料結構和函數。
pandas模組應該是python用來進行一般日常的大數據分析,最重要的模組了。
因為pandas的DataFrame資料結構,可以快速的互通於所有的二維結構資料,包括『excel,csv,json,xml,mysql,任何的資料庫,網頁的表格資料,字典dict,二維串列list』
也就是pandas的DataFrame資料結構,可以與它們之間互相簡易的存取。
然後再根據DataFrame來進行想要的大數據分析,它提供內建的演算法與資料結構,能夠用一個指令,就可以進行二維資料的『排序,篩選,關鍵字查詢,任意區間範圍查詢,統計學計算,平均值,變異數,標準差,字串檢索,字串取代,欄位樞紐分析、小記、欄位加總,把二維資料任意方式圖形視覺化顯示』
而建立pandas的DataFrame資料結構,有兩種方式,可以用column的方式來輸入資料,也可以用row的方式來輸入資料。
所以pandas是大數據分析,非常實用的利器工具,是python資料分析的首選。

(2)Numpy: Numpy專門用來處理矩陣,它的運算效率比列表list串列更高效。
Numpy是Python進行專業的數值計算之重要模組,因為大多數提供科學計算的包都是用numPy的陣列作為構建基礎,因此在進行高等數學計算時,numpy就是大數據分析的最重要工具了,因為高等數學運算,都是以矩陣的方式來進行運算,例如人工智慧,機器學習,深度學習,類神經網路計算等。

(3)sscipy:是基於numpy的科學計算包,包括統計、線性代數等工具。

(4)matplotlib:是最流行的用於繪製資料圖表的 Python 庫
也可以結合pandas模組來繪圖。

2.執行python的三種方法

1.要編寫python有三種的方法:
一、方法1:安裝python單純的python
缺點:功能陽春,沒有太多的模組,無法馬上寫大數據分析程式。
安裝網址:python官網下載

二、方法2:安裝Anaconda
優點:會同時安裝python、1000多種數學繪圖模組、Spyder編輯器,能夠支援大數據分析。
缺點:會安裝了很多你用不到的模組,浪費硬碟空間。
安裝網址:到Anacond官網下載安裝

三、方法3:使用python官網線上shell
使用repl線上python

3.Anaconda下載點

Anacond官網


3.安裝anaconda 3.安裝anaconda
功能:原始的python功能太陽春,若下載anaconda,則可以提供300多種的科學數學模組,可以提供大數據資料分析
(1)Anaconda是一個免費的Python和R語言的發行版本,用於計算科學(資料科學、機器學習、巨量資料處理和預測分析)
(2)因為Anaconda有很多的數據分析模組,所以大數據分析會使用到的『pandas、Numpy、Scipy』python package套件,在anaconda安裝完成時就已經包含在裡面了。
(3)Anaconda中文是森蚺(大蟒蛇)。
1)可以把Anaconda當作是Python的懶人包,除了Python本身(python2, 3) 還包含了Python常用的資料分析、機器學習、視覺化的套件
2).完全開源和免費
3).額外的加速、優化是收費的,但對於學術用途可以申請免費的 License
4).全平台支持:Linux、Windows、Mac
5).支持 Python 2.6、2.7、3.3、3.4,可自由切換,
6).內帶spyder 編譯器(還不錯的spyder編譯器)
7).自帶jupyter notebook 環境 (就是網頁版的python編輯器,副檔名為IPthon)

(4)常用套件:

Numpy: Python做多維陣列(矩陣)運算時的必備套件,比起Python內建的list,Numpy的array有極快的運算速度優勢
Pandas:有了Pandas可以讓Python很容易做到幾乎所有Excel的功能了,像是樞紐分析表、小記、欄位加總、篩選
Matplotlib:基本的視覺化工具,可以畫長條圖、折線圖等等…
Seaborn:另一個知名的視覺化工具,畫起來比matplotlib好看
SciKit-Learn: Python 關於機器學習的model基本上都在這個套件,像是SVM, Random Forest…
Notebook(Jupyter notebook): 一個輕量級web-base 寫Python的工具,在資料分析這個領域很熱門,雖然功能沒有比Pycharm, Spyder這些專業的IDE強大,但只要code小於500行,用Jupyter寫非常方便,Jupyter也開始慢慢支援一些Multi cursor的功能了,可以讓你一次改許多的變數名稱
(5)優點:省時:一鍵安裝完90%會用到的Python套件,剩下的再用pip install個別去安裝即可
(6)缺點:占空間:包含了一堆用不到的Python的套件(可安裝另一種miniconda)

(7)下載網址:https://www.anaconda.com/
選擇個人版:indivisual
https://www.anaconda.com/products/individual
→Download
→Windows
Python 3.7(會自動幫你安裝Python 3.7)
64-Bit Graphical Installer (466 MB)
32-Bit Graphical Installer (423 MB)

(8)安裝過程,要勾選
不勾選:add the anaconda to the system PATH(但是2020年,ananconda不建議勾選這個,容易發生錯誤)
勾選:Register anaconda as system Python 3.7

(9)安裝結束
→在windows開始→anaconda有6個項目,最常用的有3個
(1)anaconda prompt:可以直接下cmd指令
(2)Spyter:編譯器(還不錯的spyder編譯器)
(3)jupyter notebook(網頁版的python編輯器,副檔名為IPthon)
4.Anaconda prompt:cmd指令 4.使用anaconda prompt:直接下cmd指令
注意:windows 10 必須使用管理員來執行(點選anaconda prompt→滑鼠右鍵→以系統管理員身份進行)
(1)列出目前已經安裝的anaconda的模組與版本:
conda list

(2)對某個模組更新安裝
conda update 模組
範例:conda update ipython

(3)安裝某個模組
方法1:conda install 模組
範例:conda install numpy

# 安裝 NumPy 1.15 以後、 1.16 以前
conda install 'numpy>=1.15,<1.16'

方法2:pip install 模組
範例:pip install numpy

(4)解除安裝某個模組
方法1:conda uninstall 模組
範例:conda uninstall numpy

方法2:輸入 conda remove PACKAGE_NAME可以從目前的工作環境移除指定套件。
# 移除 NumPy
conda remove numpy numpy-base

方法3:pip uninstall 模組
範例:pip uninstall numpy

(5)在anaconda prompt執行python程式
方法1:
先到工作目錄:cd ch1
執行.py程式:python test1.py

方法2:python c:\chp1\test1.py

(6)常用指令
conda --version 檢視 conda 版本
conda update PACKAGE_NAME更新指定套件
conda --help 檢視 conda 指令說明文件
conda list --ENVIRONMENT 檢視指定工作環境安裝的套件清單
conda install PACAKGE_NAME=MAJOR.MINOR.PATCH 在目前的工作環境安裝指定套件
conda remove PACKAGE_NAME 在目前的工作環境移除指定套件
conda create --name ENVIRONMENT python=MAIN.MINOR.PATCH 建立新的工作環境且安裝指定 Python 版本
conda activate ENVIRONMENT 切換至指定工作環境
conda deactivate 回到 base 工作環境
conda env export --name ENVIRONMENT --file ENVIRONMENT.yml 將指定工作環境之設定匯出為 .yml 檔藉此複製且重現工作環境
conda remove --name ENVIRONMENT --all 移除指定工作環境
使用 conda list | grep numpy 檢查 Python 套件清單中是否還有 NumPy 套件
輸入 conda search PACKAGE_NAME可以檢視指定套件在 conda 中可安裝的版本列表。
# 檢視 NumPy 在 conda 中可安裝的版本
conda search numpy=1.16.3
5.用Anaconda prompt來建立虛擬環境

5.使用Anaconda prompt來建立虛擬環境
功能:可以建立多個Anaconda虛擬環境
例如:目前安裝後預設是python 3.x版本的環境,若要創建一個python 2.x的環境,就可以在Anaconda虛擬環境實現
(1)# 檢視電腦中可使用與目前所在的工作環境
conda env list

(2)使用 conda create --name ENVIRONMENT python=MAIN.MINOR.PATCH 指令可以建立出乾淨、極簡且資源隔絕的工作環境。
指令:conda create -n 虛擬環境名稱 python=版本 anaconda

# 建立一個名稱為 demo 的 Python 2 工作環境
conda create --name demo python=2
範例:建立py27env環境
conda create -n py27env python=2.7 anaconda

(3)輸入 conda activate ENVIRONMENT 可以啟動指定工作環境、
方法1:conda activate ENVIRONMENT
方法2:activate ENVIRONMENT
範例:activate py27env

方法3:到windows→開始→點選Anaconda prompt(py27env)

(4)關閉虛擬目錄,回到原本pytohn環境(base)
使用 conda deactivate 則是切換回預設的 base 工作環境。
方法1:conda deactivate
方法2:deactivate

(5)# 檢視 demo 工作環境中的套件
conda list -n py27env

(5)範例 A.建立py27env虛擬環境
conda create -n py27env python=2.7 anaconda
B.切換到py27env虛擬環境
activate py27env
C.檢視 demo 工作環境中的套件
conda list -n py27env
D.# 檢視 Python 版本
python --version
E.關閉虛擬目錄,回到原本pytohn環境(base)
deactivate

(5)複製一個與目前pyhon環境(或是py27env) 完全相同的工作環境
conda create -n 新虛擬環境名稱 --clone root
範例:conda create -n py27env2 --clone root

# 檢查明確所有虛擬環境名稱
conda info -e

(6)移除某個虛擬環境
conda remove -n 虛擬環境名稱 --all
範例:conda remove -n py27env --all

(7)常用指令整理
安裝:conda install
更新:conda update
移除:conda remove

在工作環境管理透過
創建:conda create
啟動:conda activate
停止:conda deactivate
匯出設定檔:conda env export
移除:conda remove
6.使用Spyter編譯器

6.使用Spyter:編譯器
(1)新增一個py檔案
File→ New file

print("你好,歡迎光臨")
print(1+1)

Run➤

(2)開啟已經存在的檔案
方法1:File→ Open
方法2:拖曵檔案總管的py檔案到Spyder

(3)在Spyter使用簡易智慧輸入
方法:按『tab』
範例:
先輸入p
然後按『tab』
出現list清單,都是p開始的指令

(4)程式除錯
方法1:若是這一行有指令寫錯,就會在最左邊出現三角形▲警告icon
方法2:在這個一行最左邊double click,就會出現中斷點(或是這一行按F12)
7.jupyter notebook網頁版的python編輯器

7.jupyter notebook
(1)功能:是網頁版的python編輯器,副檔名為IPthon
會開啟瀏覽器:http://localhost:8888/tree
對應的硬碟目錄 = C:\Users\電腦名稱
(例如: C:\Users\user)

(2)練習線上編輯一個簡單python程式
A.右方→New→Python3
在cell裡面輸入In[1]
a = ("apple","grape","banana")
print(a[2])

B.Run

C.修改檔案名稱→Untitled→exp1-3

D.查詢雲端檔案放置位置:C:\Users\電腦名稱\exp1-3.ipynb

(3)二種不同的Run方式
A.Run:會新增一個new cell
B.Ctrl+Enter:會停留在原本的cell

(4)在jupyter notebook使用簡易智慧輸入
方法:按『tab』
範例:
先輸入p
然後按『tab』
出現list清單,都是p開始的指令

(5)在jupyter notebook編輯的檔案無法讓python IDE編譯
jupyter notebook編輯的檔案是.ipynb
與python的.py不同
改善方法:只能把程式碼複製貼上,在兩個平台交流
8.其它線上雲端可編譯的python平台

8.其它線上雲端可編譯的python平台
網站:http://rep.it/languages/python3
 

 
chp2.python的數學函數庫:math,sympy,cmath,scipy,numpy
目錄 1.python的數學函數庫:math,sympy,cmath,scipy,numpy 範例2-1:math函數 範例2-3a:sympy數學函數 範範例2-3b:使用numpy+matplotlib, 繪圖sin(x),cos(x)
範例2-4:使用sympy繪圖sin(x) 範例2-5:使用sympy繪圖:拋物線y=x**2,10的指數函數 範例2-6:使用sympy繪圖方程式:橢圓函數x**2+y**2 = 4 範例2-7:建立符號變數(x),(x,y),(x1,x2,x3...x5) 範例2-8:建立函數表達式

1.python的數學函數庫:math,sympy,cmath,scipy,numpy

1.python的數學函數庫:math,sympy,cmath,scipy,numpy
(1)math:
功能:乘方、開方、對數,冪函數與對數函數,三角函數,角度轉換,雙曲函數,特殊函數,常量,cos,sin,e,log,tan,pow
Math模組
pi 數字常量,圓周率
e 表示一個常量
sqrt(x) 求x的平方根
fabs(x) 返回x的絕對值
factorial(x) 取x的階乘的值
fmod(x,y) 得到x/y的餘數,其值是一個浮點數
pow(x, y) 返回x的y次方,即x**y
isfinite(x) 如果x是正無窮大或負無窮大,則返回True,否則返回False
isinf(x) 如果x是正無窮大或負無窮大,則返回True,否則返回False
hypot(x) 如果x是不是無窮大的數字,則返回True,否則返回False
isnan(x) 如果x不是數字,則返回True,否則返回False
ldexp(x) 返回x*(2**i)的值////
log(x, a) 返回x的自然對數,以a為基數(不寫則預設以e為基數),a引數給定時,將x的對數返回給定的a,計算式為:log(x)/log(a
log10(x) 返回x的以10為底的對數
log1p(x) 返回x+1的自然對數(基數為e)的值
log2(x) 返回x的基2對數
modf(x) 返回由x的小數部分和整數部分組成的元組
trunc(x) 返回x的整數部分
ceil(x) 取大於等於x的最小的整數值,如果x是一個整數,則返回x
floor(x) 取小於等於x的最大的整數值,如果x是一個整數,則返回自身
radians(x) 把角度x轉換成弧度,與degrees 為反函式
degrees(x) 把x從弧度轉換成角度
sin(x) 求x(x為弧度)的正弦值
sinh(x) 求x(x為弧度)的雙曲正弦值
cos(x) 求x的餘弦,x必須是弧度
tan(x) 返回x(x為弧度)的正切值
tanh(x) 返回x(x為弧度)的雙曲正切值
copysign(x, y)把y的正負號加到x前面,可以使用0
exp(x) 返回math.e,也就是2.71828的x次方
expm1(x) 返回math.e的x(其值為2.71828)次方的值減1
frexp(x) 返回一個元組(m,e),其計算方式為:x分別除0.5和1,得到一個值的範圍
fsum(x) 對迭代器裡的每個元素進行求和操作
gcd(x,y) 返回x和y的最大公約數

(2)sympy:
SymPy是一個符號計算的Python庫。它的目標是成為一個全功能的計算機代數系統,同時保持代碼簡潔、易於理解和擴展。它完全由Python寫成,不依賴於外部庫。[2]
SymPy支持符號計算、高精度計算、模式匹配、繪圖、解方程、微積分、組合數學、離散數學、幾何學、概率與統計、物理學等方面的功能

python有強大的科學計算庫,包括Numpy,Scipy,Ipython,matplotlib,pandas,sympy等等。
其中只有sympy是可以使用代數及其他各種數學符號來作運算的。
(其他庫偏向於各種資料值運算和工程計算),
因為sympy的開發目的就是成為一個全功能的代數系統,這就比較合純數學的計算(高中數學,大學數學)

Sympy是一個數學符號庫(sym代表了symbol,符號),包括了積分,微分方程,三角等各種數學運算方法,是工科最基本的數學函數庫,用起來媲美matlab,而且其精度比math函數庫精確。
例如:
多項式的合併、展開、化簡
微積分運算
微分方程求解
線性代數運算
係數匹配
繪圖
simplify運算式化簡,solve方程自動求解,limit求極限,diff求導,dsolve()計算微分方程
intergrate積分計算:1.定積分,2.不定積分,3.雙重定積分,4. 雙重不定積分

進階:SymPy 可以支援的內容包括但不限於:
1. 範疇論(Category Theory);
2. 微分幾何(Differential Geometry);
3. 常微分方程(ODE);
4. 偏微分方程(PDE);
5. 傅立葉轉換(Fourier Transform);
6. 集合論(Set Theory);
7. 邏輯計算(Logic Theory)。

(3)cmath:專門用來處理複數運算。

(4)scipy:
SciPy是一個開源的Python演算法庫和數學工具包。
SciPy包含的模組有最佳化、線性代數、積分、插值、特殊函數、快速傅立葉變換、訊號處理和圖像處理、常微分方程式求解和其他科學與工程中常用的計算。
與其功能相類似的軟體還有MATLAB、GNU Octave和Scilab。

(5)numpy:
NumPy是Python語言的一個擴充程式庫。支援高階大量的維度陣列與矩陣運算,此外也針對陣列運算提供大量的數學函式函式庫。
NumPy的核心功能是ndarray(即n-dimensional array,多維陣列)資料結構。這是一個表示多維度、同質並且固定大小的陣列物件。而由一個與此陣列相關聯的資料型態物件來描述其陣列元素的資料格式
NumPy提供了與MATLAB相似的功能與操作方式,因為兩者皆為直譯語言,並且都可以讓使用者在針對陣列或矩陣運算時提供較純量運算更快的效能
NumPy也可以結合其它的Python擴充函式庫。例如SciPy,這個函式庫提供了更多與MATLAB相似的功能;以及Matplotlib,這是一個與MATLAB內建繪圖功能類似的函式庫。

Numpy和Scipy是機器學習專案中最受歡迎的函式庫
在GitHub上的ML專案中,7成以上的機器學習專案都使用了Numpy函式庫,Scipy函式庫則是占了近5成
Numpy函式庫具備平行處理能力,支援多維度陣列和矩陣的數學運算,在機器學習專案中,能夠處理龐大的資料量,因此以超過7成的機器學習專案占有率,成為最受歡迎的ML函式庫,
另外,負責運算的Scipy、管理資料集Pandas和提供視覺化的Matplotlib函式庫,在ML專案中,都有超過4成的占有率。

範例2-1:math函數

2.範例2-1:math函數
import math

#印出所有math函數名稱
print(dir(math))

#常數pi, e
print('常數 pi=', math.pi)
print('常數 e=', math.e)

#三角函數
#角度 = 30度
#degrees(x) 傳回弧度 x 的角度 (degree)
#radians(x) 傳回角度 x 的弧度 (radian)
d1 = 30
r1 = math.radians(d1)
print('sin(30度) = ', math.sin(r1))

#反三角函數: arc Sin(0.5) = pi/6
r2 = math.asin(0.5)
d2 = math.degrees(r2)
print('ashin(0.5)= radians =',r2)
print('ashin(0.5)= degrees =',d2)

#factorial(x) 傳回 x 階乘 =x!
d3 = math.factorial(5)
print('5!=', d3)

#絕對值
d4 = abs(-9)
d5 = math.fabs(-7)
print('絕對值 = abs(-9)=', d4)
print('絕對值 = math.fabs(-7)=', d5)

# 自然指數 = math.exp(1)
print('自然指數 = math.exp(1) = ', math.exp(1))
print('自然指數 = math.e = ', math.e)

#求20,16的最大公因數
print('求20,16的最大公因數=', math.gcd(20,16))

#四捨五入,無條件捨去,無條件進位
print('6.7的四捨五入 =', round(6.7))
print('6.7的無條件捨去 =', math.floor(6.7))
print('6.7的無條件進位 =', math.ceil(6.7))

#5的3次方
print('5的3次方=', 5**3)
print('5的3次方=', math.pow(5,3))

#開根號
print('25開根號 = ', math.sqrt(25))

成果圖示:
程式碼內容

範例2-3a:sympy數學函數

3.範例2-3a:sympy數學函數
import sympy,math

#印出所有sympy函數名稱
print(dir(sympy))

#sympy比math的精度高:sympy的精度比math高
print('math,5開根號**2 = ', math.sqrt(5)**2)
print('sympy,5開根號**2 = ',sympy.sqrt(5)**2)

#常數pi, e
print('常數 pi=', sympy.pi)
print('sin(pi/6) = ', sympy.sin(sympy.pi/6))
#圓周率(pi)與尤拉常數(e),在 SymPy 中以符號的方式進行運算,若是要算出實際的浮點數,可以使用 evalf() 這個函數
print('sym.pi的實際數字=sympy.pi.evalf()=', sympy.pi.evalf())
print('常數 E=sympy.E.evalf()=', sympy.E.evalf())

#三角函數
#角度 = 30度
r1 = math.radians(30)
d1 = math.degrees(sympy.pi/6)
print(r1,d1)
#d1 = sympy.degree((sympy.pi)/6)
print('math, sin(30度) = ', math.sin(sympy.pi/6))
print('sympy, sin(30度) = ', sympy.sin(sympy.pi/6))
print('sympy, sin(30度) = ', sympy.sin(r1))

#反三角函數: arc Sin(0.5) = pi/6
r2 = sympy.asin(0.5)
d2 = math.degrees(r2)
print('ashin(0.5)= radians =',r2)
print('ashin(0.5)= degrees =',d2)

#factorial(x) 傳回 x 階乘 =x!
d3 = sympy.factorial(5)
print('5!=', d3)

#絕對值
d5 = sympy.Abs(-7)
print('絕對值 = sympy.Abs(-7)=', d5)

# 自然指數 = sympy.exp(1)
print('自然指數 = sympy.exp(1) = ', sympy.exp(1))

#求20,16的最大公因數
print('求20,16的最大公因數=', sympy.gcd(20,16))

#四捨五入,無條件捨去,無條件進位ceiling
print('6.7的四捨五入 =', round(6.7))
print('6.7的無條件捨去 =', sympy.floor(6.7))
print('6.7的無條件進位 =', sympy.ceiling(6.7))

#5的3次方 sympy.Pow()
print('5的3次方=', 5**3)
print('5的3次方=', sympy.Pow(5,3))

#開根號
print('25開根號 = ', sympy.sqrt(25))

成果圖示:
程式碼內容

範範例2-3b:使用numpy+matplotlib, 繪圖sin(x),cos(x)

4.範例2-3:使用numpy+matplotlib, 繪圖sin(x),cos(x)
import numpy as np
import matplotlib.pyplot as plt

#一個圖:sin(x)
#numpy的pi = np.pi
x = np.arange(0,np.pi*8,0.01)
#numpy 的sin = np.sin
y = np.sin(x)
plt.plot(x,y,color='red')
plt.show()

#二個圖:sin(x), cos(x)
x = np.arange(-np.pi*4,np.pi*4,0.01)
y1 = np.sin(x)
y2 = np.cos(x)
plt.plot(x,y1,color='red',label='y=sinx')
plt.plot(x,y2,color='blue',label='y=cosx')
plt.title = 'sin(x) and cos(x)'
plt.xlabel('angel(pi)')
plt.ylabel('amplitude',rotation=0)
#圖例標籤
plt.legend()
#plt.legend(['y=sinx','y=cosx'])
plt.show()

#三個圖:sin(x), cos(x), sin(x)+cos(x)
x = np.arange(-np.pi*4,np.pi*4,0.01)
y1 = np.sin(x)
y2 = np.cos(x)
y3 = y1 + y2
plt.plot(x,y1,color='red',label='y=sinx')
plt.plot(x,y2,color='blue',label='y=cosx')
plt.plot(x,y3,color='green',label='y=sinx+cosx')
plt.title = 'sin(x) and cos(x)'
#設定xlabe,ylabell的位置:position=
plt.xlabel('angel(pi)',position=(0.95,0))
plt.ylabel('amplitude',position=(0,1),rotation=0)

#圖例標籤的位置 loc = "lower left"
'''
'best' 0
'upper right' 1
'upper left' 2
'lower left' 3
'lower right' 4
'right' 5
'center left' 6
'center right' 7
'lower center' 8
'upper center' 9
'center' 10
'''
plt.legend(loc="lower left")
#plt.legend(['y=sinx','y=cosx'])

#設定:不顯示圖片的四邊黑框線
ax = plt.gca()
ax.spines["top"].set_color("none")
ax.spines["right"].set_color("none")
ax.spines["bottom"].set_position(("data", 0))
ax.spines["left"].set_position(("data", 0))

plt.show()

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

範例2-4:使用sympy繪圖sin(x)

5.範例2-4:使用sympy繪圖sin(x)
from sympy import *
#設定變數x
x = symbols('x')
plot(sin(x))

#指定x範圍畫圖
plot(sin(x),(x,-5, 5))

#畫三條圖:y=x, y=sin(x), y=cos(x)
plot(x,sin(x),cos(x),(x,-5, 5))

#設定畫圖顏色
p1 = plot(x,sin(x),cos(x),(x,-5, 5),show=False)
p1[0].line_color = 'green'
p1[1].line_color = 'r'
p1.show()

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

範例2-5:使用sympy繪圖:拋物線y=x**2,10的指數函數y=10**x,標準指數函數y=e**x

6.範例2-5:使用sympy繪圖:拋物線y=x**2,10的指數函數y=10**x,標準指數函數y=e**x
from sympy import *
x = symbols('x')

#拋物線:y=x**2
plot(x**2,(x,-2, 2))

#10的指數函數:y=10**x
plot(10**x,(x,-2, 2))

#標準指數函數:y=e**x
plot(E**x,(x,-2, 2))

#三個圖畫在一起
p1 = plot(x**2,10**x,E**x,(x,-2, 2),show=False)
#設定顏色
p1[0].line_color = 'r'
p1[1].line_color = 'g'
#設定圖例legend
p1[0].label = '$x^2$'
p1[1].label = '$10^x$'
p1[2].label = '$e^x$'
#要打開 True
p1.legend = True

p1.show()

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

範例2-6:使用sympy繪圖方程式:橢圓函數x**2+y**2 = 4,雙曲線函數:x**2/4-y**2 =1

7.範例2-6:使用sympy繪圖方程式:橢圓函數x**2+y**2 = 4,雙曲線函數:x**2/4-y**2 =1
from sympy import *
#隱函數:implicit function
#x**2+y**2 = 4,可以轉換成 y=sqrt(4-x**2)
#若方程式,有雙變數,就要設定雙變數:x, y = symbols('x y')
x, y = symbols('x y')
#建立:橢圓函數:x**2+y**2 = 4
plot_implicit(Eq(x**2 + y**2,4),(x,-3,3),(y,-3,3))

#建立:雙曲線函數:x**2/4-y**2 =1
plot_implicit(Eq(x**2/4-y**2, 1))

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

範例2-7:建立符號變數(x),(x,y),(x1,x2,x3...x5)

8.範例2-7:建立符號變數(x),(x,y),(x1,x2,x3...x5)
from sympy import *

#建立符號變數x
x = symbols('x')
print('x=',x)

#建立符號變數x,y(其中 x >0)
x,y = symbols('x y')
x = symbols('x', positive = True)
print('x=',x)

#建立5個變數,x1,x2,x3,x4,x5
x = symbols('x1:6')
print('x1,....x5=',x)

成果圖示:
程式碼內容

範例2-8:建立函數表達式

8.範例2-8:建立函數表達式
from sympy import *
txt = 'x**2 + 2*x + 5'
formular = sympify(txt)
print('formular=', formular)
成果圖示:
程式碼內容

 

 
chp3.多項式,方程式,函數
1.多項式,方程式,函數的定義與差異 2.如何在python的sympy印出標準的數學符號 3.以下的範例,最好在jupyter編輯器上執行 範例3-1:多項式= (x+2y+3)**2 範例3-2:多項式的因式分解,因式拆解
範例3-3:方程式f(x)=0,求解x=? (solve,roots) 7.方程式,函數,畫圖的方法不一樣:plot(函數),plot_implicit(Eq方程式) 範例3-4:二元一次方程式聯立求解 範例3-5:兩個多項式聯立求解(圓,與線的交叉點 9.sympy的各種繪圖指令
範例3-6:3D曲面函數(z=f(x,y),不是方程式,是函數) 範例3-7:三元一次方程式f(x,y,z)聯立求解(圓,與線的交叉點)

1.多項式,方程式,函數的定義與差異:

1.多項式,方程式,函數的定義與差異:
(1).多項式,polynomial
定義:由文字符號(通常為x)和數字,進行加法和乘法運算所構成的代數式。
例如:x**2+3x+1
例外:若代數式的文字符號(x)出現在分母、根號裡,或絕對值中,則此代數式就不能定義為「多項式」。
(例如:1/x,sqrt(x),則不能稱為多項式)

(2).方程式,Equation
定義:方程式乃是令(多項式 = 0)的式子,
特色:方程式,可以解方程式,可以求出滿足此式子的 x 的值為多少

(3).函數,Function
定義:只要在集合x中的每一個元素,都能在f集合,找到剛好一個元素與之對應,就形成函數對應
例如:f(x) = x**2+3x+1
特色:只要給定一個x值,就可以算出其對應的f

(4).sympy如何處理『多項式,方程式,函數』
■多項式,polynomial
指令:poly(x**2+2*x+1)

■方程式,Equation
指令:Eq(x**2+2*x+1, 0)

■函數,Function
指令:
x = symbols('x')
f = Function('f')(x)
f = x**2+2*x+1
求值:
num = f.subs(x,2)

2.如何在python的sympy印出標準的數學符號

2.如何在python的sympy印出標準的數學符號
方法:
#先import可以輸出在畫面顯示方程式的latex格式
from sympy.interactive import printing
printing.init_printing(use_latex=True)
from sympy import *

#再display數學符號
f = poly(x**2+2*x+1)
display(f)

#注意:若要顯示標準數學符號最好用jupyter
若在spider編輯器,display則會顯示在plot視窗

#若是用print()輸出,則是純文字格式的數學符號

3.以下的範例,最好在jupyter編輯器上執行

3.以下的範例,最好在jupyter編輯器上執行
jupyter的優點:
(1)在Jupyter不需要print()指令,按 Run就能把執行結果顯示出來

(2)在Jupyter若要顯示標準數學符號,也不需要
from sympy.interactive import printing
printing.init_printing(use_latex=True)
display(f)

而是直接在in[]裡面單獨輸入:f,如何按Run,就能夠顯示標準數學符號

(3)在jupyter裡面弱若是忘記某個指令的參數用法語法,可以
直接在in[]裡面單獨輸入:expr.series?
就會顯示語法說明

範例3-1:多項式= (x+2y+3)**2

4.範例3-1:多項式= (x+2y+3)**2
#可以輸出在畫面顯示方程式的latex格式
from sympy.interactive import printing
printing.init_printing(use_latex=True)
from sympy import *

x = symbols('x')

#一元多項式
f = poly(x**2+2*x+1)
display(f)

#一元多項式,相乘,展開
f = poly((x+3)*(x+5))
print(f)

#多項式相除
g = poly((x**3+x**2+x+1)/(x-2))
display(g)

#二元多項式
x,y = symbols('x y')
f = poly('(x+2*y+3)**2')
display(f)

#一元函數求對應值:f(x=2)
#求值f(x=3),兩種寫法:(1) f.subs({x:3}。(2)f.subs(x,3)
f = Function('y')(x)
f = (x+2)**2
print('f(x=2)=',f.subs({x:3}))
print('f(x=2)=',f.subs(x, 3))

#二元函數求值:g(x=2,y=1)
g = Function('g')(x,y)
g = (x+2*y+3)**2
print('g(x=2,y=1) = ', g.subs({x:2, y:1}))


成果圖示:
程式碼內容

範例3-2:多項式的因式分解,因式拆解

5.範例3-2:多項式的因式分解,因式拆解
from sympy import *
x = symbols('x')
f1= poly(x**2-7*x+10)

#多項式的:因式分解:factor(f1,x)
#factor(expr, x) 函數可以處理因式分解
print('多項式的因式分解=', factor(f1, x))
print('多項式的因式分解=', factor(x**2-7*x+10, x))

#因式拆解 = apart(f)
f2 = poly(1/((1 + x)* (3 + x)))
f3 = apart(f2)
print('多項式的因式拆解=', apart(f2))

#因式合併 = together(p)
print('多項式的因式合併=', together(f3))

#若要顯示標準方程式,用display()
from sympy.interactive import printing
printing.init_printing(use_latex=True)
display(f2)

成果圖示:
程式碼內容

範例3-3:方程式f(x)=0,求解x=? (solve,roots)

6.範例3-3:方程式f(x)=0,求解x=? (solve,roots)
from sympy import *
x = symbols('x')
f1 = Eq(x**2-7*x+10, 0)
#方程式f1(x)=0,求解x=?

#方法1:roots(f)求解
print('f1(x)=0,求解x=?',roots(f1))

#方法2:solve(f)求解
print('f1(x)=0,求解x=?',solve(f1))
print('f1(x)=0,求解x1=',solve(f1)[0])
print('f1(x)=0,求解x2=',solve(f1)[1])

#畫出f的圖形
#若要畫出圖形,必須是函數f(x),才能畫圖
f2 = Function('f2')(x)
f2 = x**2-7*x+10
plot(f2,(x,-3,8))

#方程式f2(x)=x**2-9*x+10=0,求解x=?
f2 = Eq(x**2-9*x+10,0)
print('f2(x)=0,求解x=?',solve(f2))
print('f2(x)=0,求解x1=',solve(f2)[0])
print('f2(x)=0,求解x2=',solve(f2)[1])

成果圖示:
成果圖示:

標準寫法程式碼
極簡化寫法程式碼

7.方程式,函數,畫圖的方法不一樣:plot(函數),plot_implicit(Eq方程式)

7.方程式,函數,畫圖的方法不一樣:plot(函數),plot_implicit(Eq方程式)
(1)方程式畫圖:plot_implicit(Eq方程式)
#函數f(x)=...的畫圖用:plot(f1,line_color='red')
#方程式的畫圖用:plot_implicit(eq1,line_color='red')
#語法:plot_implicit(expr, x_var=None, y_var=None, adaptive=True, depth=0, points=300, line_color='blue', show=True, **kwargs
#預設只能畫出一個圖
#錯誤指令:plot_implicit(f1, f2)

(2)進階:同時畫出兩個方程式equation的方法
#使用plot_implicit(Eq()畫出兩個圖的方法:使用extend(p2),最後再p1.show())
p1 = None
p1 = plot_implicit(f1,show=False)
p2 = plot_implicit(f2,show=False)
p1.extend(p2)
p1.show()

(2)函數畫圖的方法:plot(函數)
f2 = Function('f2')(x)
f2 = x**2-7*x+10
plot(f2,(x,-3,8))

範例3-4:二元一次方程式聯立求解

範例3-4:二元一次方程式聯立求解
from sympy import *
x, y = symbols('x y')
#二元一次聯立方程式
f1 = Eq(x + 2*y - 8, 0)
f2 = Eq(2*x - y - 6, 0)

#方程式畫出圖
#函數f(x)=...的畫圖用:plot(f1,line_color='red')
#方程式的畫圖用:plot_implicit(eq1,line_color='red')
#語法:plot_implicit(expr, x_var=None, y_var=None, adaptive=True, depth=0, points=300, line_color='blue', show=True, **kwargs
#預設只能畫出一個圖
#錯誤指令:plot_implicit(f1, f2)

#進階:同時畫出兩個方程式equation的方法
#使用plot_implicit(Eq()畫出兩個圖的方法:使用extend(p2),最後再p1.show())
p1 = None
p1 = plot_implicit(f1,show=False)
p2 = plot_implicit(f2,show=False)
p1.extend(p2)
p1.show()

#求解兩條線的交叉點(solve(f1,f2))
print('解二元一次聯立方程式=', solve((f1,f2),(x,y)))

成果圖示:
成果圖示:
標準寫法程式碼
極簡化寫法程式碼

範例3-5:兩個多項式聯立求解(圓,與線的交叉點)

8.範例3-5:兩個多項式聯立求解(圓,與線的交叉點)
from sympy import *
x, y = symbols('x y')
#二元一次聯立方程式
f1 = Eq(x**2 + y**2-12, 0)
f2 = Eq(2*x - y - 6, 0)

#畫出圖
#使用plot_implicit(Eq()畫出兩個圖的方法:使用extend(p2),最後再p1.show())
p1 = None
p1 = plot_implicit(f1,(x,-5,5),(y,-5,5),show=False)
p2 = plot_implicit(f2,(x,-5,5),(y,-5,5),show=False)
p1.extend(p2)
p1.show()

#求解兩條線的交叉點(solve(f1,f2))
print('解二元一次聯立方程式=', solve((f1,f2),(x,y)))

from sympy.interactive import printing
printing.init_printing(use_latex=True)
display(solve((f1,f2),(x,y)))

成果圖示:
成果圖示:
標準寫法程式碼
極簡化寫法程式碼

9.sympy的各種繪圖指令

9.sympy的各種繪圖指令
(1)各種繪圖指令:
plot: Plots 2D line plots.
plot_parametric: Plots 2D parametric plots.
plot_implicit: Plots 2D implicit and region plots.
plot3d: Plots 3D plots of functions in two variables.
plot3d_parametric_line: Plots 3D line plots, defined by a parameter.
plot3d_parametric_surface: Plots 3D parametric surface plots.

(2)限制:若要畫出3D圖,目前只能畫出2個變數的3D圖(plot3d)
10.sympy的繪圖功能,目前的限制與適用狀況如下:
plot: Plots 2D line plots.
plot_parametric: Plots 2D parametric plots.
plot_implicit: Plots 2D implicit and region plots.
plot3d: Plots 3D plots of functions in two variables.
plot3d_parametric_line: Plots 3D line plots, defined by a parameter.
plot3d_parametric_surface: Plots 3D parametric surface plots.

所以,無法畫出3個變數的3D圖

(3)二次曲面的種類:


範例3-6:3D曲面函數(z=f(x,y),不是方程式,是函數)

10.範例3-6:3D曲面函數(z=f(x,y),不是方程式,是函數)
from sympy import *
#from sympy.plotting import plot3d

#(1).橢圓拋物面 z = x^2/a^2 + y^2/b^2
#錯誤方法:用方程式Eq的觀念,用plot_implicit來畫3D圖
#結果:失敗,因為plot_implicit只能畫出2D圖
x,y,z = symbols('x y z')
f = Eq(x**2/2**2 + y**2/2**2 -z, 0)
#plot_implicit(f)

#修正方法,使用 z=f(x,y)的函數方式,來畫圖
x,y = symbols('x y')
z = Function('z')(x,y)

z = x**2/2**2 + y**2/2**2
plotting.plot3d(z)
print('橢圓拋物面求值(x=5,y=5)=', z.subs({x:5,y:5}))

#(2).雙曲拋物面 z = x^2/a^2 - y^2/b^2
x,y = symbols('x y')
z = Function('z')(x,y)
z = x**2/3**2 - y**2/4**2
plotting.plot3d(z)
print('雙曲拋物面求值(x=5,y=5)=', z.subs({x:5,y:5}))

#單葉雙曲面 f = x^2/a^2 + y^2/b^2 -z^2/c^2- 1
x,y,z = symbols('x y z')
f = Function('f')(x,y,z)
f = x**2/3**2 + y**2/2**2 -z**2/2**2 -1
#這個函數,無法畫出3D圖
#因為plot3d()最多只能畫出兩個變數的圖形,這個題目的變數有三個
#plotting.plot3d(f)
#定義三個symbols,可以計算函數值
print('單葉雙曲面求值(x=5,y=5,z=5)=', f.subs({x:5,y:5,z:5}))

成果圖示:
成果圖示:
成果圖示:
成果圖示:
標準寫法程式碼
極簡化寫法程式碼

範例3-7:三元一次方程式f(x,y,z)聯立求解(圓,與線的交叉點)

11.範例3-7:三元一次方程式f(x,y,z)聯立求解(圓,與線的交叉點)
from sympy import *
x,y,z = symbols('x y z')
f1 = Eq(2*x - y + z - 10, 0)
f2 = Eq(3*x + 2*y - z - 16, 0)
f3 = Eq(x + 6*y - z - 28, 0)

#注意:solve求解的格式,有兩種:
#(1)solve([f1,f2,f3])。其中要注意[]
#(2)solve((f1,f2,f3),(x,y,z))
#print('三個多項式f(x,y,z)聯立求解=', solve([f1,f2,f3]))
print('三個多項式f(x,y,z)聯立求解=', solve((f1,f2,f3),(x,y,z)))

成果圖示:
標準寫法程式碼
極簡化寫法程式碼

 
 
chp4.複數complex number,與極坐標系統(Polar coordinate system)
1.極坐標系統(Polar coordinate system) 範例4-1:極坐標系統(Polar coordinate system) 範例4-2:複數系統cmath的應用:三種坐標轉換(複數坐標,極坐標,直角坐標) 範例4-3:虛數的物理意義, 範例4-4:在極坐標上畫出一個點,畫出一條線
範例4-5:畫出極坐標的曲線方程式 範例4-6:用sympy + plot3d_parametric_line畫3D曲線/a> 範例4-7:用sympy + plot3d_parametric_line畫3D曲面(surface)
1.極坐標系統(Polar coordinate system)

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))

成果圖示:

程式碼內容

 
 
chp5.微分,函數的微分
目錄 1.一階微分,一階導數 範例5-1:求微分,常見的幾類微積分函數基本公式 範例5-2:求在某點位置的微分值(斜率) 4.二階微分,二階導數
範例5-3:求二階導數,判別反曲點,判別線型凹凸 範例5-4:求一階導數,判別局部極大值,極小值的位置

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')

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

 
chp6.積分,integral
1.積分定義 2.比較:積分,函數,微分,二階微分的關係 3.sympy的積分語法:integrate 範例6-1:求函數的不定積分,定積分
範例6-2:比較sympy與Maxima軟體在計算積分的準確率 範例6-3:三角函數的積分

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

成果圖示:
程式碼內容


 
 
chp7.微分方程式,ODE:ordinary differential equation
1.常微分方程式(ordinary differential equation,簡稱ODE) 2.求解一階微分方程式的唯一解,必須給定一個邊界條件 3.求解二階微分方程式的唯一解,必須給定兩個邊界條件 4.sympy的微分函數寫法有兩種
範例7-1:微分方程式ODE求解,y'= ∂y/∂x = x 範例7-2:微分方程式ODE求解,y'(x) = x**3+1 範例7-3:微分方程式ODE求解,y'(x)-2x/y=0 範例7-4:微分方程式ODE求解,y'(x) = 3x**2(y+2)
範例7-5:試求在 xy 平面中通過點(1,1)且斜率為-y/x之曲線

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')

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

 
 
chp8.二階常係數微分方程式,Second-Order Differential Equations
範例8-1:二階常係數齊性 ODE 的解法(一)--相異實根 範例8-2:二階常係數齊性ODE的解法(二)--重根 範例8-3:二階常係數齊性ODE的解法(三)--複數根 範例8-4:Euler-Cauchy 方程式之解--複數根
範例8-5:二階常係數非齊性ODE之特解 範例8-6:高階(三階)常係數齊性 ODE 之求解 範例8-7:一階聯立微分方程式,求特解 範例8-8:二階聯立微分方程式,求通解(y1,y2)

範例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)

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

 
 
chp9.矩陣算法,解聯立方程式,線性代數
1.聯立方程式的矩陣求法 範例9-1:聯立方程式的聯立求法(線性代數)

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()

成果圖示:
程式碼內容

 
 
chp10.泰勒級數,泰勒多項式,用泰勒展開式模擬函數,Taylor Series
1.泰勒展開:以多項式來逼近函數 2.泰勒級數的應用-1:大約估值/a> 3.泰勒一次項係數即導數 4.用微分計算泰勒多項式的係數
5.多項式函數的圖形特色 6.典型範例:由銷售量q與單價x的三次多項式關係,求月營收的極大值極小值 7.完整的n次泰勒多項式函數(以導數表示係數),及其發生極值的位置 範例10-1:使用sympy把多項式轉成泰勒級數
範例10-2:f(1.97)泰勒級數估計值大概為-76.03,請計算精確值 範例10-3:求三次多項式的極值發生地方 範例10-4:由銷售量q與單價x的三次多項式關係,求月營收的極大值極小值 2.正餘弦sin(x),cos(x)的泰勒多項式
13.指數函數的泰勒多項式 13.對數函數的泰勒多項式 範例10-5:以泰勒多項式逼近指數函數exp(x) 範例10-6:以泰勒多項式逼近sin(x)函數
範例10-7:以泰勒多項式逼近ln(x)函數

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:讀取網址的網頁


 
 
chp11.傅立葉級數,Fourier Series
1.傅立葉級數的應用 2.傅立葉數學的學習內容與步驟 3.傅立葉級數公式 4.傅立葉級數的物理意義
5.傅立葉級數公式的視覺化Series 範例11-1:在(-pi < x < pi)內用傅立葉級數模擬三次方多項式的線型 範例11-2:比較不同n級數的傅立葉級數去模擬三次方多項式的 範例11-3:比較不同n級數的傅立葉級數去模擬指數函數exp(x)
範例11-4:比較不同n級數的傅立葉級數去模擬自然對數函數ln(x) 範例11-5:求週期函數Square Wave方波函數f(x)的傅立葉級數F 範例11-6:求週期函數f(x)=x**2,在(-pi < x < pi)的傅立葉級數 範例11-7:求週期2L函數f(x)=|x|,在(-3 < x <3)的傅立葉級數
範例11-8:求週期Square Wave方波函數f(x)的傅立葉級數 範例11-9:傅立葉積分:求週期為∞的方波函數f(x),)的傅立葉積分

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

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

 
 
chp12.傅立葉轉換(Fourier Transform, FT)
1.何謂傅立葉變換 2.傅立葉轉換公式有兩種 3.傅立葉變換:就是兩個坐標的變換 4.Inverse Fourier Transform,反向傅立葉轉換
5.什麼時候要用到傅立葉變換? 6.傅立葉轉換的可視化 7.傅立葉級數,傅立葉轉換的關聯與差異 8.傅立葉變換家族
9.傅立葉級數(FS) vs 傅立葉轉換(FT)的應用差異 10.傅立葉積分 vs 傅立葉轉換的差異 11.波的表示式有兩種:歐拉公式:e^(ix) = cos(x) + i*sin(x) 12.傅立葉轉換(FT)後是個複數: g(ω) = A(ω)+ iP(ω)
13.傅立葉轉換的應用:相機的美顏功能,與變聲 14.python有三種模組都支援快速傅立葉轉換 範例12-1:求週期為∞的指數函數f(x)=exp(-x**2)的傅立葉轉換 範例12-2:D.T.F.T:離散資料,非週期數據,的傅立葉轉換法FFT
範例12-3:D.T.F.T:使用spicy轉換離散資料,非週期數據的傅立葉轉換法FFT 範例12-4:把雜訊濾波:D.T.F.T:使用scipy做離散非週期性數據的傅立葉轉換 範例12-5:從2個sin合成波訊號中濾除雜訊

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()

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

檔案


 
 
chp13.網絡大數據爬取與分析5(Selenium自動化網頁操作)
1.日期時間函式庫 範例13-1:模組datetime(顯示日期,時間 範例13-2:模組date(只能顯示日期) 範例13-3:模組timedelta(間隔日期,時間) 範例13-4:傳回自1970年1月1日凌晨0:0:0開始至今的秒數
範例13-5:計算呼叫某個函式所需要的時間 2.可以迭代(iterate)的函式庫= itertools,函式enumerate,zip,filter,map,reduce 範例13-6:顯示陣列 範例13-7:zip(a,b)是結合兩個可迭代物 範例13-8:模組itertools:給定起始值,遞增值,就可以產生無窮數列
範例13-9:filter(過濾函數def或lambda,資料列) 範例13-10:過濾模組map 範例13-11:操作每個元素模組reduce 範例13-12:套件collections-模組OderedDict:可將兩個串列,合併形成成一個字典dict 範例13-13:模組deque:操作可迭代資料列,新增,刪除,反轉
例13-14:模組dCounter:可計算可迭代資料列相同元素的次

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())
程式碼內容


 
 
chp14.數學函數庫math,sympy,微分,積分,偏微分
1.存取xml:模組:xml.etree.ElementTree 範例14-1:讀取person.xml的所有每個節點資訊,查詢所有的mail,查詢卓水信資料 範例14-2:讀取person.xml的所有每個節點資訊,查詢所有的mail,查詢卓水信 範例14-3:修改並存入xml
線上XML/JSON互相轉換工具 2.讀取網頁:request(url) 範例14-4:讀取網頁:web = request.urlopen(網址) 3.存取 json(模組:json)
範例14-5:轉成json:jumps。轉成dict:loads 範例14-6:讀取網絡上的json檔案 範例14-7:讀取電腦上的json檔案 3.網路爬蟲BeautifulSoup
範例14-8:讀取網頁標題 範例14-9:讀取網址的網頁

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)
程式碼內容