用 Python+SciPy+SymPy+GeoGebra 執行微積分計算 1

賴博伯

已于 2024-03-28 11:53:25 修改

阅读量1.6k

点赞数

分类专栏：上課講義笔记筆記文章标签： python 矩阵算法 GeoGebra 几何学

于 2021-05-26 17:22:25 首次发布

本文链接：https://blog.csdn.net/m0_47985483/article/details/117298121

版权

笔记同时被 3 个专栏收录

47 篇文章 0 订阅

订阅专栏

上課講義

46 篇文章 1 订阅

订阅专栏

筆記

32 篇文章 1 订阅

订阅专栏

“Talk is cheap. Show me the code.”
― Linus Torvalds

老子第41章
上德若谷
大白若辱
大方無隅
大器晚成
大音希聲
大象無形
道隱無名

拳打千遍, 身法自然

GeoGebra 在中學數理科到大二微積分相當夠用
Python 中執行矩陣運算的主力是 NumPy
Python 中執行數值計算的主力是 SciPy(先載入NumPy)
Python 中執行微積分或符號運算的主力是 SymPy
符號計算相當於叫電腦模擬我們中學或大學時用手算出數學公式的計算

如果考量到

免費開源社群穩定
兼具數值與符號運算
又能延伸到研究所

建議使用 GeoGebra 同時搭配 Python+(先載入)NumPy(矩陣)+SciPy(科學計算)+SymPy(符號運算,
(如果是GeoGebra 搭配以下之軟體,可能數值與符號運算就無法兼具,

GeoGebra 搭配以下之免費軟體, : GeoGebra + R (偏重統計與數值), GeoGebra + Octave (偏重數值), GeoGebra + Maxima (偏重符號運算) 等,
如果不在意付費, 則搭配付費軟體, 學習曲線會較低, 因為付費軟體通常說明文檔會較齊全貼心: GeoGebra + Matlab (偏重數值), GeoGebra + Maple (偏重符號運算), GeoGebra + Mathematica (偏重符號運算) 搭配, 都可.)

建議初學者使用何種軟體教微積分或計算微積分

最好是同時使用 GeoGebra 與一個程式語言類型的軟體

例如 GeoGebra 搭配 Python(+Scipy(NumPy)+SymPy)
依照本人教授微積分多年, 嘗試用各種軟體製作教材同時並教會同學使用, 以輔助同學自我學習, 我目前得到的心得是:
最好是同時使用 GeoGebra 與一個程式語言類型的軟體,
例如 GeoGebra 搭配 Python(+Scipy(先載入NumPy)+SymPy),
或 GeoGebra 搭配 Matlab, Maple, Mathematica 或 Octave, R, Maxima 等搭配, 都可,
Matlab, Maple, Mathematica 是付費的, 另外近年, 還有Julia 網路常介紹.

我自己的建議是可以兩方搭配一起用, 當然學習的時間, 以我本人的經驗, 教師以課餘瑣碎的時間邊教邊摸索自學, 約需1~3年的時間(單只有GeoGebra, 在有熟悉的人帶領學習或在旁邊可以隨時發問的前提下只需2,3 個月, 如果完全自己摸索沒人可以問, 可能要半年才能製作較進階的GeoGebra 內容, 至於 Python+NumPy+SciPy+SymPy 等, Python 雖然好學, 但是要掌握這麼多的工具箱, 還是要不短的時間, 如果完全自己摸索沒人可以問, 假設連 Python 也是從頭學, 可能至少要半年, 才能製作基礎的教材, 閒暇有空就學, 應該 1 年多就可以進入熟悉的狀態, 加上前面的 GeoGebra 半年的時間, 合計約至少兩年差不多, 另外 Python+NumPy+SciPy+SymPy 如果有人教或有跟隨課程應該半年就可以基本夠用了), 但是, 作者觀察本校師範院校教出的學生, 尤其是我自己教過的同學進入各中小學教書, 如果不肯把握時間加緊自學(中小學教學上處理學生的瑣事很多, 時不時有一堆表面功夫的研習必須參加, 應校長要求申請計畫, 應付上級督導等, 必須是對使用軟體輔助數學探索有熱愛且學習意志力很堅定的老師), 過了10年, 發現這個學生還是不會駕馭軟體協助自己教學與研究, 只能羨慕自己的同事, 輕易嫻熟的操控軟體做出精美的數學教材, 電腦這樣的東西, 花時間敲指令下去就有收穫, 所謂 "拳打千遍, 身法自然 ", 如果用看的用想的, 永遠只能羨慕別人.

GeoGebra 在中學數理科到大一微積分相當夠用

GeoGebra 可以直接以滑鼠拖拉畫圖, 以滑桿功能展示即時動態變化(動畫), 也有不少微積分的指令, 例如以下我用 GeoGebra 製作的積分的動態展示, 加上文字說明, 按鈕等, 同樣的效果, 如果僅是靜態圖, 用 Python+Matplotlib 的 pyplot畫, 恐怕要花上兩三倍的時間.
而且下圖之 Geogebra 的圖是可以拖動滑桿, 改變切割區間之粗細, 動態演示的, 這種動畫效果, 如果要用 Python+Matplotlib 的 pyplot 畫呈現, 會花上更多倍時間, 且需要對 Matplotlib 有更多的了解與熟悉.
(更正: 下圖中的說明 “長方形法” 其實是"梯形法")
數值積分梯形法講解用 Sequence
例如以下是我用 GeoGebra 製作的微積分考卷的圖, 要呈現 Varberg 微積分那本原文教科書的某題之圖, 加上自己的設計, 用 GeoGebra 製作 3D 函數圖也很輕鬆:
201606光通期末考極值1
或, 以下是用 GeoGebra 製作的兩變數函數 $f(x,y)=4+xy-x^2-y^2$ 侷限在單位圓上求極值, 的說明圖:
Varberg_P670_12.9Example5_20210413jpg
可以按右鍵, 轉動視角,
同樣的效果, 如果僅是靜態圖, 用 pyplot, 恐怕要手忙腳亂一陣子, 3D 可能還要動用到 MayaView 等, 如果是動畫, 還得查一下 Tkinter 的指令等.

另一方面, 如果完全用 GeoGebra 算微積分, 因為 GeGebra 主力是動態幾何, 會覺得指令不夠, 要想有完整的微積分, 最佳化, 矩陣, 解微分方程等函數, 還是要搭配依賴較完整的程式語言, 例如 Python+SciPy(先載入NumPy)+SymPy, 或 Matlab, Maple, Mathematica , 或 R, Maxima 等.

所以我自己的建議是可以兩方搭配一起用, 當然學習的時間, 以我本人的經驗, 教師以課餘瑣碎的時間邊教邊摸索自學, 約需1~3年的時間(有熟悉的人帶領學習的前提下, GeoGebra只需2,3個月, 自己摸索可能要半年才能製作較進階的GeoGebra內容, 另外 Python+NumPy+SciPy+SymPy 如果有人教或有跟隨課程應該半年就可以基本夠用了).

微積分軟體之建議搭配	是否免費	動態幾何(滑鼠+微積分指令+3D)	程式語言	數值與符號運算功能
	皆是	GeoGebra	Python+SciPy(先載入NumPy)+SymPy	皆有
	皆是	GeoGebra	R	數值與統計運算較強
	皆是	GeoGebra	GNU Octave	數值運算較強(指令類似Matlab)
	皆是	GeoGebra	Maxima	符號運算較強
	一半	GeoGebra	Matlab	數值運算較強
	一半	GeoGebra	Maple	符號運算較強
	一半	GeoGebra	Mathematica	符號運算較強
註	皆是	desmos	desmos 與 GeoGebra 的繪圖計算器非常類似, 雖然簡易上手, 但是不像GeoGebra功能較全面性

大家可以看的出來, 初學者乾脆選第一項: GeoGebra 搭配 Python(+SciPy(先載入NumPy)+SymPy), 既免費開源, 又數值計算與符號運算功能一次到位,
不過學習過程較長, 不像商用軟體, 所有指令都一次到位, Python(+SciPy(先載入NumPy)+SymPy) 需多花時間熟悉, 那些指令是 Python 內建的, 那些在 scipy, 那些在 sympy, 算是一個小的困擾.
作者用過 Maple 多年, Maple是以符號運算知名, 深知, 符號運算較強的軟體在數值運算遠不如Matlab 之指令豐富與效能快速,
作者也用過 Matlab 多年, , Matlab 是以矩陣數值運算知名, 深知, 數值運算較強的軟體在符號運算遠不如 Maple之指令豐富,
而偏偏, 微積分計算的特質, 就是混合數值計算與符號運算才能完整,
以作者教學的經驗, 在大一微積分範圍, 可以只用符號運算強的軟體就很夠用, 甚至只有 GeoGebra也很夠了, 故如果初學者選 GeoGebra 搭配符號運算強的軟體 Maxima, Maple, Mathematica 等, 或數值計算強的軟體 Matlab, Octave 等都可以.
但是如果考量後續的延伸性, 延伸到數學系大二大三的數值計算, 科學計算, 數學建模, 資料分析等, 工程學系的高等工程數學, 運籌學系的單純形法, 動態規劃法, 非線性規劃等等, 選擇搭配符號運算強的軟體會開始感到效能不濟, 尤其在跑一些尺度稍大的演算法時, 就會走上作者的同樣的路, 先學 Maple, 之後必須再學 Matlab 才夠用.
故如果考量後續的延伸大二以上到研究所, 建議初學者選擇 GeoGebra 搭配數值計算強的軟體,
當然如果一開始就選 GeoGebra 搭配 Python(+SciPy(先載入NumPy)+SymPy), 既免費開源, 又數值計算與符號運算功能一次到位, 就不用煩惱這些.

是否免費	大一大二微積分範圍	延伸至大二以上的科學計算,非線性規劃,高等工程數學等	數值與符號運算功能
皆是	GeoGebra+ Python+SciPy(先載入NumPy)+SymPy	GeoGebra+ Python+SciPy(先載入NumPy)+SymPy	皆有
Octave, C, Java 是	GeoGebra+數值計算強的軟體 Matlab, Octave R, C, Java 等	GeoGebra+數值計算強的軟體 Matlab, Octave, R, C, Java 等	數值運算較強
Maxima是	GeoGebra+符號運算強的軟體 Maxima, Maple, Mathematica 等	大尺度數值模擬較耗時	符號運算較強
是	GeoGebra		數值與符號運算功能皆夠大一大二用

desmos 與 GeoGebra 的比較

近幾年網路似乎開始鼓吹 desmos 一款類似 GeoGebra 的軟體, 較簡易上手, 但是我發現功能相當受限, 介面與 GeoGebra 的繪圖計算器類似, 不知是誰學誰的?
兩個都可以畫 3D 函數圖.

跟 GeoGebra 比起來, desmos 的功能只是 GeGebra的一小部分, 但是 desmos 似乎較會做宣傳.

以下是 GeoGebra官網線上的繪圖計算器, 可以把他想成是 GeoGebra 6 的一個特例的介面狀況
GeoGebra的繪圖計算器

以下是 desmos 官網線上的介面, 與 GeoGebra 官網的線上繪圖計算器非常類似:

desmos 官網線上的介面

Python 中有數學指令的模組或程式庫

以下大略介紹 Python 中有數學指令的模組或程式庫, 部分參考自陳擎文教學網, https://acupun.site/lecture/python_math/index.htm#chp5 link

內建模組或外部程式庫	包含內容	大一大二微積分範圍	延伸至大二以上的科學計算,非線性規劃,高等工程數學等	數值與符號運算功能
內建(原生) math 模組	基本數學函數與指令: pi, e, sqrt(x) , fabs(x), factorial(x), fmod(x,y), sin(x), sinh(x), cos(x), tan(x)等	不足應付	不足應付	數值運算
內建(原生 cmath	處理複數運算	不足應付	不足應付	數值運算
矩陣運算 NumPy 外部程式庫	NumPy 支援高維度, 大尺度矩陣(即array)運算，此外也針對矩陣運算提供大量的數學函式函式庫。在矩陣計算範圍內, NumPy 提供了與 MATLAB 相似的功能與操作方式	不適合 (較適合線性代數課程)	支援高維度矩陣運算	數值運算
數值運算 SciPy (底層含 NumPy) 外部程式庫	SciPy 包含的模組有最佳化、線性代數、積分、插值、特殊函數、快速傅立葉變換、訊號處理和圖像處理、常微分方程式求解和其他科學與工程中常用的計算。與其功能相類似的軟體還有MATLAB、GNU Octave和Scilab。	足夠應付微積分數值計算部分	足夠應付數值運算部分	數值運算
符號運算 SymPy 外部程式庫	Sympy是一個數學符號庫（sym代表了symbol，符號），包括了積分，微分方程，三角等各種數學運算方法，是工科最基本的數學函數庫，用起來媲美Matlab, Maple，而且其精度比math函數庫精確。例如：多項式的合併、展開、化簡微積分運算微分方程求解線性代數運算係數匹配繪圖進階：SymPy 可以支援的內容包括但不限於： 1. 範疇論（Category Theory）； 2. 微分幾何（Differential Geometry）； 3. 常微分方程（ODE）； 4. 偏微分方程（PDE）； 5. 傅立葉轉換（Fourier Transform）； 6. 集合論（Set Theory）； 7. 邏輯計算（Logic Theory）。	足夠應付	足夠應付所有符號運算與小尺度數值運算部分	符號運算與部分數值運算

Python 中執行微積分的符號運算的主力是 SymPy

特別指出, 經過作者在網路搜尋, 發現,
一般在網路講解用 Python 算微積分的網誌或課程, 主要還是用 SymPy 來做各種微積分的計算,
SymPy 是提供求符號計算的程式庫, 類似 Maple, Mathematica 等, 但是不像Maple, Mathematica是需付費的, SymPy是開源免費,只需在 Python 裝好之後再加裝第三方程式庫進來就好, 如果是使用 Anaconda 或 Google的 CoLab, 則都會已經自動裝好. (Matlab 裡搭配符號計算的引擎是 MuPad, 最一開始是用 Maple引擎, 不知何原因後來改成用Mupad).

Python 可以執行 calculus 的程式庫:

scipy: 可以提供求~~數值的微分~~, 數值的積分, 零根, min, max, 解 ODE 等數值計算, scipy 已包含 NumPy 在內.
sympy: 可以提供求符號計算的微分, 積分, Taylor series, 零根, min, max, 解 ODE 等符號計算, 當然數值計算也能, 在單單微積分範圍內, sympy的指令還是較齊全!, 但是學習符號計算通常較困難, 需要較多的數學素養.

對於學生, 如果只需要計算數值解的微分, 積分, 求零根, 極大極小值(min, max), 解常微分方程ODE等等, 可以只學 SciPy 即可, SymPy 的學習曲線較陡峭.

數值計算與符號計算之差別

所謂　數值微分 numerical differentiation、數值積分 numerical integration,
是求其近似值, 而非公式.
例如以下 cos(x) 的數值積分是求其近似值, 而非公式:
$\int_0^{1.2} cos(x) dx = 0.9320390859672262$

>>> from scipy import integrate

>>> integrate.quad(scipy.cos,0,1.2)
(0.9320390859672262, 1.0347712530914001e-14)

相對於此是符號計算 symbolic computation 的微分或積分, 也就是我們大學用手算的微分或積分的解析解, 有明確的數學式子的解答,
(注意, 不是所有的積分都有解析解, 例如橢圓積分, 或是例如 $\int sin(x^2) dx$ 等, 都沒有明確的解析式子的表法, 對於這類的積分, 只能求近似解(數值解), 或是用特殊函數表達其解),

符號計算是求公式, 例如以下 cos(x) 的手算積分是求公式而非近似值:
$\int cos(x) dx = sin(x) +C$

也就是符號計算相當於叫電腦模擬我們中學或大學時用手算出數學公式的計算

以下是 SymPy 的語法

>>> import sympy
>>> x = sympy.Symbol('x')
>>> x
x
>>> sympy.integrate( sympy.cos(x), x)
sin(x)
>>> sympy.diff( sympy.cos(x), x)
-sin(x)

註: math 模組與 numpy 都沒數值微分、數值積分的指令, numpy.diff() 是求 array差值.
Python 中似乎缺乏數值微分的指令, 通常都是用 sympy 來算微分,
scipy 似乎也只有數值積分,
讀者可以自己寫一個求數值微分的函數, codes 可以參考後面數值微分那節.

以下陳擎文教授此站主要介紹使用 sympy 算微積分的課程講義:
陳擎文教學網：python求解數學式（高中數學，大學數學，工程數學，微積分）
方程式，多項式，函數，三角函數，微分，積分，泰勒級數，傅立葉級數，傅立葉積分，傅立葉轉換，一階微分方程式，二階微分方程式，三階微分方程式，偏微分方程式，聯立方程式矩陣解法，張量，向量的坐標轉換，複數，極坐標
https://acupun.site/lecture/python_math/index.htm#chp5 link

Python, NumPy, SciPy, SymPy 的安裝或線上使用

安裝Python
請參考本人的另一篇
https://blog.csdn.net/m0_47985483/article/details/109522800
安裝Python 那節
link
安裝 NumPy, SciPy, Matplotlib, Pandas
Python安裝之後並沒有 NumPy, SciPy, Matplotlib, Pandas, 他們是額外加裝在 Python 上的程式庫,
在 Windows 下, 打開 “命令提示字元” 的視窗, 輸入

>> pip install numpy
>> pip install scipy
>> pip install sympy

或是使用 Anaconda, 安裝好之後, 最重要的程式庫都已裝好,
Anaconda + Jupyter Notebook 會自動安裝好所需的科學計算或大數據的程式庫 (or Anaconda + Spyder or Anaconda + PyCharm 等),
線上使用可以用 Google Colab, 也會自動安裝好所需的科學計算或大數據的程式庫.

Python 中執行數值計算的主力是 SciPy (需先載入NumPy)

SciPy 的底層是 NumPy.
scipy 的入門說明講義不太好找, 常常會搜尋到別的頁面,
在 https://scipy-lectures.org/index.html這裡,
link

SciPy的模組種類

https://docs.scipy.org/doc/scipy/reference/index.html

scipy的模組很多, 本篇只關注跟微積分較相關的兩個
Integration and ODEs (scipy.integrate)
Optimization and root finding (scipy.optimize)

API Reference
The exact API of all functions and classes, as given by the docstrings. The API documents expected types and allowed features for all functions, and all parameters available for the algorithms.

Clustering package (scipy.cluster)
Constants (scipy.constants)
Discrete Fourier transforms (scipy.fft)
Legacy discrete Fourier transforms (scipy.fftpack)
Integration and ODEs (scipy.integrate)
Interpolation (scipy.interpolate)
Input and output (scipy.io)
Linear algebra (scipy.linalg)
Miscellaneous routines (scipy.misc)
Multidimensional image processing (scipy.ndimage)
Orthogonal distance regression (scipy.odr)
Optimization and root finding (scipy.optimize)
Signal processing (scipy.signal)
Sparse matrices (scipy.sparse)
Sparse linear algebra (scipy.sparse.linalg)
Compressed sparse graph routines (scipy.sparse.csgraph)
Spatial algorithms and data structures (scipy.spatial)
Special functions (scipy.special)
Statistical functions (scipy.stats)
Statistical functions for masked arrays (scipy.stats.mstats)
Low-level callback functions

sciPy 求零根

1.6.5.3. Finding the roots of a scalar function

ref: https://scipy-lectures.org/intro/scipy.html#finding-the-roots-of-a-scalar-function link

用 scipy.optimize.root(fun, ,,)求零根

scipy.optimize.root(fun, x0, args=(), method=‘hybr’, jac=None, tol=None, callback=None, options=None)

官網的例子
$f(x) = x^2 + 10\sin(x)$
求其零根.

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt

def f(x):
    return x**2 + 10*np.sin(x)

x = np.arange(-10, 10, 0.1)

plt.plot(x, f(x))

root = optimize.root(f, x0=1)
print(root)
print('root1=', root.x)

root2 = optimize.root(f, x0=-2.5)
print('root2=', root2.x)

輸出:

>>> 
= RESTART: C:\Users\user\Desktop\1.6.5.3. Finding the roots of a scalar function.py
    fjac: array([[-1.]])
     fun: array([0.])
 message: 'The solution converged.'
    nfev: 10
     qtf: array([1.33310463e-32])
       r: array([-10.])
  status: 1
 success: True
       x: array([0.])
root1= [0.]
root2= [-2.47948183]

SciPy 求數值微分、數值積分

所謂　數值微分 numerical differentiation、數值積分 numerical integration, 是求其近似值,
相對於此, 是符號計算 symbolic computation 的微分或積分,
也就是我們大學用手算的微分或積分的解析解, 有明確的數學式子的解答.

數值微分

math 模組與 numpy 都沒數值微分、數值積分的指令, numpy.diff 是求 array差值.
Python 模組與程式庫中似乎缺乏數值微分的指令, 通常都是用 sympy 來算微分,
scipy 似乎只有數值積分,
數值微分的原理很簡單, 照定義算, 取極限的意思, 就是切夠小格就好, 可以自己寫一個求數值微分的函數.

Ex: 自己寫一個求數值微分的函數.
Sol: 自己寫一個求數值微分的函數並不需要用到 scipy, 只需用 Python 或至多 numpy,
根據微分的定義:
$f'(x_0) = \lim_{ \nabla x \rightarrow 0 } \frac{f(x_0+\nabla x) - f(x_0)}{\nabla x}$
故只要用夠小的 $\nabla x$ 去計算
$\frac{f(x_0+\nabla x) - f(x_0)}{\nabla x}$
就是 $f'(x_0)$ 的近似值.

以下自定義函數 `numericalDiff(f, x_0, n)`

輸入引數中的 f, 為求微分之函數, 可以用 lambda 表達式 輸入,
x_0 為求微分之點,
n, 代表切 1/n 這麼小, i.e., $\nabla x = 1/n$ ,
如果要精度較高, 可以輸入較大的 n, 例如 10000, 代表 $\nabla x = \frac{1}{10000}=0.0001$ .

lambda 表達式 :
例如
lambda x: x**2
代表一個臨時產生的函數 $f(x) = x^2$ .

import math
lambda x: math.sin(x**2)
(注意: sin(x) 需要使用內建的 math 模組)
代表一個臨時產生的函數 $f(x) = sin(x^2)$ .

>>> import math
>>> f = lambda x: math.sin(x**2)
>>> f(2)
-0.7568024953079282

自定義數值微分的函數 numericalDiff() 如下方 codes,
呼叫此函數之用法, 示例:
numericalDiff( lambda x: x**2, 1, 100)
代表對函數 $f(x) = x^2$ 求微分, 在 x=1 的位置, $\nabla x$ 取為 $\frac{1}{100}=0.01$ .
如果要精度較高, 可以輸入較大的 n,
例如 n=10000, 代表 $\nabla x = \frac{1}{10000}=0.0001$ ,

自定義數值微分的函數 numericalDiff() 的程式碼:

# By Prof. P-J Lai MATH NKNU 20210605
# Numerical differentiation
# numericalDiff.py

def numericalDiff(f, x_0, n):
    deltaX = 1/n
    #print( 'the Df of x^2 at {} is {}'.format(x_0, (f(x_0+deltaX) - f(x_0))/deltaX)  )
    print( f'the Df of x^2 at {x_0} is { (f(x_0+deltaX) - f(x_0))/deltaX }' )

    return (f(x_0+deltaX) - f(x_0))/deltaX

numericalDiff( lambda x: x**2, 1, 100)

輸出:

>>> 
================ RESTART: C:/Users/user/Desktop/numericalDiff.py ===============
the Df of x^2 at 1 is 2.0100000000000007

數值積分

當我們想要用 scipy 算數值積分時,
需先下 from scipy import integrate, 載入 scipy 之 integrate 模組,
(或使用 from scipy import * 載入所有指令)
不能只載入 scipy, 再直接使用 scipy.integrate.quad(), 例如以下

>>> import scipy
>>> scipy.integrate.quad(scipy.cos,0,1.2)
Traceback (most recent call last):
  File "<pyshell#1>", line 1, in <module>
    scipy.integrate.quad(scipy.cos,0,1.2)
AttributeError: module 'scipy' has no attribute 'integrate'

以下才正確
先下 from scipy import integrate
再呼叫 integrate.quad(scipy.cos,0,1.2),

>>> from scipy import integrate

>>> integrate.quad(scipy.cos,0,1.2)
(0.9320390859672262, 1.0347712530914001e-14)

或 scipy.integrate.quad(scipy.cos,0,1.2) 皆可

>>> from scipy import integrate
>>> scipy.integrate.quad(scipy.cos,0,1.2)
(0.9320390859672262, 1.0347712530914001e-14)

查詢 scipy.integrate 中的積分指令

>>> from scipy import integrate
>>> dir(scipy.integrate)
Traceback (most recent call last):
  File "<pyshell#10>", line 1, in <module>
    dir(scipy.integrate)
NameError: name 'scipy' is not defined
>>> dir(integrate)
['AccuracyWarning', 'BDF', 'DOP853', 'DenseOutput', 'IntegrationWarning', 'LSODA', 'OdeSolution', 'OdeSolver', 
'RK23', 'RK45', 'Radau', 
'__all__', '__builtins__', '__cached__', 
'__doc__', '__file__', '__loader__', '__name__', '__package__', '__path__', '__spec__', '_bvp', '_dop', '_ivp', 
'_ode', '_odepack', '_quad_vec', '_quadpack', '_quadrature', 'complex_ode', 'cumtrapz', 'cumulative_trapezoid', 
'dblquad', 'fixed_quad', 'lsoda', 'newton_cotes', 'nquad', 'ode', 'odeint', 'odepack', 'quad', 'quad_explain', 
'quad_vec', 'quadpack', 'quadrature', 'romb', 'romberg', 'simps', 'simpson', 'solve_bvp', 'solve_ivp', 
'test', 'tplquad', 'trapezoid', 'trapz', 'vode']

官網 1.6.7.1. Function integrals 之例子

求

$\int_0^{\frac{\pi}{2}} sin(x) dx \sim 1$

res, err = quad(np.sin, 0, np.pi/2)
res 是輸出之結果, err 是誤差,
np.allclose(res, 1) 是詢問 res 與理論值 1 是否接近?
(註: math.isclose(a, b) 功能類似)

>>> from scipy.integrate import quad
>>> import numpy as np
>>> res, err = quad(np.sin, 0, np.pi/2)
>>>> res
0.9999999999999999
>>> err
1.1102230246251564e-14
>>> np.allclose(res, 1)
True
>>> np.allclose(err, 1 - res)
True
>>>

1.6.5 Optimization and fit 求極值與曲線擬合(內外插法)

Ref: Scipy Lecture Notes, 1.6.5 Optimization and fit: scipy.optimize, link

params, params_covariance = optimize.curve_fit(test_func, x_data, y_data, p0=[2,2])

#!/usr/bin/env python
# coding: utf-8

# # 6.5 Optimization and fit: scipy.optimize
# Ref: Scipy Lecture Notes — Scipy lecture notes主頁2020

# In[3]:


from scipy import optimize
import numpy as np


# In[4]:


x_data = np.linspace(-5,5, num=50)


# In[5]:


x_data


# In[6]:


y_data = 2.9*np.sin(1.5*x_data) + np.random.normal(size=50)


# In[7]:


y_data


# In[8]:


def test_func(x, a, b):
    return a*np.sin(b*x)


# In[9]:


params, params_covariance = optimize.curve_fit(test_func, x_data, y_data, p0=[2,2])


# In[10]:


print(params)
##output
##>>> 
##= RESTART: C:\data\NEW\網路免費軟體資料\Python教學\Python科學計算\
## sci_optimize_求##min_max_zeros_求零根可以轉成求min\
## 6.5 Optimization and fit scipy.optimize.py
##[2.97083835 1.52207958]


# In[11]:


import matplotlib.pyplot as plt


# In[24]:


plt.figure(figsize=(10,4))
plt.scatter(x_data, y_data, label='data')
x=np.linspace(-5,5, num=50)
y=test_func(x,params[0],params[1])
plt.plot(x,y,color='red', label='optimize.curve_fit')
plt.legend()
plt.show()


# In[ ]:

Scipy Lecture Notes, 6.5 Optimization and fit_scipy.optimize_1

(1.6.5) 2.7 scipy.optimize 求極值

scipy 的 optimizer種類:
scipy.optimize 基本上是求極植(max, min),
至於求零根, 不管事純量函數 scalar function 或向量值函數 vector-valued function
都可以將求零根問題
$F(x)=(f_1(x), f_2(x),\cdots,f_n(x))=0,$
轉為求min 的問題
$\min \left(\sum_{i=1}^n f_i(x)^2 \right).$

以下是scipy.optimize 中各種可以使用的求極植的方法(函數), 參考自
Ref: Scipy Lecture Notes, 1.6.5 Optimization and fit: scipy.optimize, link
2.7. Mathematical optimization: finding minima of functions, link
20210603 online.

這裡有簡單總結何種情況使用何種 optimizer 求解器:

2.7.5. Practical guide to optimization with scipy
2.7.5.1. Choosing a method 選擇適合的 optimizer 求解器 link

result = optimize.minimize_scalar(f)
optimize.minimize(f, [2, -1], method="CG")
optimize.minimize(f, [2,-1], method="Newton-CG", jac=jacobian)
result = optimize.minimize(g, [1, 1], method="Powell", tol=1e-10)

optimize.minimize(f, np.array([0, 0]), method="SLSQP",
constraints={"fun": constraint, "type": "ineq"})

result = optimize.minimize_scalar(f, bracket=(-5, 2.9, 4.5), method="Brent",
options={"maxiter": iter}, tol=np.finfo(1.).eps)

optimize.minimize(f, np.array([0, 0]), method="L-BFGS-B",
bounds=((-1.5, 1.5), (-1.5, 1.5)))

optimize.minimize(f, [2, -1], method="Nelder-Mead")

optimize.brute(f, ((-1, 2), (-1, 2)))

optimize.minimize(f, np.array([0, 0]), bounds=((-1.5, 1.5), (-1.5, 1.5)))

optimize.minimize(f, x0, constraints={"fun": constraint, "type": "ineq"})

2.7.5.1. Choosing a method
All methods are exposed as the method argument of scipy.optimize.minimize().

Without knowledge of the gradient
• In general, prefer BFGS or L-BFGS, even if you have to approximate numerically
gradients. These are also the default if you omit the parameter method - depending
if the problem has constraints or bounds
• On well-conditioned problems, Powell and Nelder-Mead, both gradient-free
methods, work well in high dimension, but they collapse for ill-conditioned problems.
With knowledge of the gradient
• BFGS or L-BFGS.
• Computational overhead of BFGS is larger than that L-BFGS, itself larger than
that of conjugate gradient. On the other side, BFGS usually needs less function
evaluations than CG. Thus conjugate gradient method is better than BFGS at
optimizing computationally cheap functions.
With the Hessian
• If you can compute the Hessian, prefer the Newton method (Newton-CG or TCG).
If you have noisy measurements
• Use Nelder-Mead or Powell.

scipy 的 constants 模組

ref: http://yhhuang1966.blogspot.com/2020/04/python-scipy.html link

scipy 的 constants 必須先 import

只有 sci.pi, sci.e 等是 scipy 內建的,
其餘常數都必須先
from scipy import constants
才能叫出 constants下的常數 golden, c等

>>>import scipy
>>> scipy.pi
3.141592653589793
>>> scipy.c
Traceback (most recent call last):
  File "<pyshell#10>", line 1, in <module>
    scipy.c
AttributeError: module 'scipy' has no attribute 'c'
>>> scipy.golden
Traceback (most recent call last):
  File "<pyshell#11>", line 1, in <module>
    scipy.golden
AttributeError: module 'scipy' has no attribute 'golden'
>>> scipy.e
2.718281828459045

如果已載入scipy, 想直接用 scipy.constants.golden 是不行的
必須特別單獨載入 constants,
from scipy import constants
才看的到 constants 的屬性, 也才能用 dir(constants)

以下是錯誤的

>>>import scipy
>>> scipy.constants.golden
Traceback (most recent call last):
  File "<pyshell#1>", line 1, in <module>
    scipy.constants.golden
AttributeError: module 'scipy' has no attribute 'constants'

>>> dir(scipy.constants)
Traceback (most recent call last):
  File "<pyshell#2>", line 1, in <module>
    dir(scipy.constants)
AttributeError: module 'scipy' has no attribute 'constants'
>>> dir(constants)
Traceback (most recent call last):
  File "<pyshell#3>", line 1, in <module>
    dir(constants)
NameError: name 'constants' is not defined

以下是正確的

>>> from scipy import constants
>>> dir(constants)
['Avogadro', 'Boltzmann', 'Btu', 'Btu_IT', 'Btu_th', 'ConstantWarning', 'G', 'Julian_year', 'N_A', 'Planck', 'R', 'Rydberg', 'Stefan_Boltzmann', 
'Wien', '__all__', '__builtins__', '__cached__', '__doc__', '__file__', '__loader__', '__name__', '__package__', '__path__', '__spec__', 
'_obsolete_constants', 'acre', 'alpha', 'angstrom', 'arcmin', 'arcminute', 'arcsec', 'arcsecond', 'astronomical_unit', 'atm', 'atmosphere',
 'atomic_mass', 'atto', 'au', 'bar', 'barrel', 'bbl', 'blob', 'c', 'calorie', 'calorie_IT', 'calorie_th', 'carat', 'centi', 'codata', 'constants', 
 'convert_temperature', 'day', 'deci', 'degree', 'degree_Fahrenheit', 'deka', 'dyn', 'dyne', 'e', 'eV', 'electron_mass', 'electron_volt', 
 'elementary_charge', 'epsilon_0', 
 'erg', 'exa', 'exbi', 'femto', 'fermi', 'find', 'fine_structure', 'fluid_ounce', 'fluid_ounce_US', 
 'fluid_ounce_imp', 'foot', 'g', 'gallon', 'gallon_US', 'gallon_imp', 'gas_constant', 'gibi', 'giga', 'golden', 'golden_ratio', 
 'grain', 'gram', 'gravitational_constant', 
 'h', 'hbar', 'hectare', 'hecto', 'horsepower', 'hour', 'hp', 'inch', 'k', 'kgf', 'kibi', 'kilo', 'kilogram_force', 
 'kmh', 'knot', 'lambda2nu', 'lb', 'lbf', 'light_year', 'liter', 'litre', 'long_ton', 'm_e', 'm_n', 'm_p', 'm_u', 'mach', 'mebi', 'mega', 
 'metric_ton', 'micro', 'micron', 'mil', 'mile', 'milli', 'minute', 'mmHg', 'mph', 'mu_0', 'nano', 'nautical_mile', 'neutron_mass', 
 'nu2lambda', 'ounce', 'oz', 'parsec', 'pebi', 'peta', 'physical_constants', 'pi', 'pico', 'point', 'pound', 'pound_force', 'precision', 
 'proton_mass', 'psi', 'pt', 'short_ton', 'sigma', 'slinch', 'slug', 'speed_of_light', 'speed_of_sound', 'stone', 'survey_foot', 
 'survey_mile', 'tebi', 'tera', 'test', 'ton_TNT', 'torr', 'troy_ounce', 'troy_pound', 'u', 'unit', 'value', 'week', 'yard', 'year', 
 'yobi', 'yotta', 'zebi', 'zepto', 'zero_Celsius', 'zetta']

>>> constants.c
299792458.0
>>> constants.golden
1.618033988749895

SymPy

從 SymPy 的模組可以看出 SymPy 主要還是在數學的深度與廣度的設計, 主力還是在數學範圍, 底下有較深的數學主題, 例如李群 Lie algebra, 範疇論Category Theory, 密碼學Cryptography, 偏微分方程 PDE, 張量 tensor 等等, 用這來教微積分, 應該是綽綽有餘, 但是對於高效能的數值算法, 她並不擅長, 還是要用 SciPy 或 Matlab 等輔助.

Ref: SymPy 的模組可以在此查看: SymPy Modules Reference, https://docs.sympy.org/latest/modules/index.html link

SymPy 的模組很多

Algebras
Introduction
Quaternion Reference
Assumptions
Calculus
Category Theory
Combinatorics
Cryptography
Diophantine
Discrete
Geometry
Inequality Solvers
Lie Algebra
Logic
Matrix Expressions
Number Theory
ODE
PDE
Physics
Plotting
Sets
Simplify
Solvers
Stats
Tensor
Vector
等等

官網操作手冊

官網手冊目前 2021 的 1.8 版 , 跟比 2018 的 1.3 版比, 就多了 7百多頁, 更新的速度相當驚人!
Ref: Welcome to SymPy’s documentation!, https://docs.sympy.org/latest/index.html link

2.2.1 What is Symbolic Computation 求 $\sqrt{8}$

Symbolic computation deals with the computation of mathematical objects symbolically. This means that the mathematical objects are represented exactly, not approximately, and mathematical expressions with unevaluated variables are left in symbolic form.

官網文件先舉開更號的例子,
在 Python 中我們要算 $\sqrt{9}$ , 算出來的是浮點數

>>> 9**(1/2)
3.0

或是用 math.sqrt()

>>> import math
>>> math.sqrt(9)
3.0

此時可以使用符號運算的功能, 來保持手寫數學式子的呈現,

>>> sympy.sqrt(5)
sqrt(5)

如果覺得 sqrt(5) 不像手寫的風格,
可以使用 sympy.pprint()美化列印, 印出為一如我們平時手寫數學式子的風格

>>> sympy.pprint(sympy.sqrt(5))
√5

此時 sympy不會立即傳回浮點值, 而是保持在一個數學符號 $\sqrt{5}$ 的’狀態
但是碰到再對他平方, 還是能夠算出為 5, 而不是一個類似字串的死符號,

>>> sympy.sqrt(5)**2
5

另一種好處是, 例如碰到 $\sqrt{8}$ , 他會回我們手寫計算的習慣呈現方式 $2\sqrt{2}$

>>> sympy.sqrt(8)
2*sqrt(2)
>>> sympy.pprint(sympy.sqrt(8))
2⋅√2

另外一種情況是, 如果我們想呈現
$(\sqrt{5})^2$
這樣的寫法, 又希望他同時能對其做相對的數學運算
如果用以上的指令, 會立刻傳回浮點值

>>> (math.sqrt(5))**2
5.000000000000001

可以使用延遲求值的方式,
(Prevent expression evaluation, Ref: Advanced Expression Manipulation, https://docs.sympy.org/latest/tutorial/manipulation.html#prevent-expression-evaluation link)
這種方式在 Maple 裡或是一般的符號運算軟體裡也常用

一種方法是可以用 sympy.UnevaluatedExpr(sympy.sqrt(5))
將他包裹起來

>>> expSqrt5 = sympy.UnevaluatedExpr(sympy.sqrt(5))
>>> expSqrt5
sqrt(5)

之後再對他平方, 也不會展開為 5:

>>> expSqrt5**2
(sqrt(5))**2
>>> sympy.pprint(expSqrt5**2)
  2
√5

如果要求值, 可以先用 doit(), 解除延遲求值

>>> expSqrt5**2.doit()
SyntaxError: invalid syntax
>>> (expSqrt5**2).doit()
5

再用 evalf()

>>> expSqrt5.doit().evalf()
2.23606797749979

另一種延遲求值方法是, sympy.sqrt(5,evaluate=False)

>>> sympy.sqrt(5,evaluate=False)
sqrt(5)

但是再對他平方, 會展開為 5

>> expSqrt5_1=sympy.sqrt(5,evaluate=False)
>>> expSqrt5_1
sqrt(5)
>>> expSqrt5_1**2
5

import sympy

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

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

x,y

(x, y)

type(x)

sympy.core.symbol.Symbol

foo=x+2*y

foo

$\displaystyle x + 2 y$

type(foo)

sympy.core.add.Add

foo+1

$\displaystyle x + 2 y + 1$

foo-y

$\displaystyle x + y$

x*foo

$\displaystyle x \left(x + 2 y\right)$

(x*foo).expand()

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

expandFoox = sympy.expand(x*foo)

expandFoox

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

sympy.factor(expandFoox)

$\displaystyle x \left(x + 2 y\right)$

2.2.3 The Power of Symbolic Computation

expr = sympy.sin(x**2)*sympy.exp(2*x)

expr

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

DiffExpr = sympy.diff(expr)
DiffExpr

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

sympy.integrate(DiffExpr,x)

$\displaystyle 2 \left(\int e^{2 x} \sin{\left(x^{2} \right)}\, dx + \int x e^{2 x} \cos{\left(x^{2} \right)}\, dx\right)$

expr1 = sympy.sin(x)*sympy.exp(x)
expr1

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

DiffExpr1 = sympy.diff(expr1,x)

DiffExpr1

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

sympy.integrate(DiffExpr1)

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

sympy.integrate(sympy.sin(x**2),x)

$\displaystyle \frac{3 \sqrt{2} \sqrt{\pi} S\left(\frac{\sqrt{2} x}{\sqrt{\pi}}\right) \Gamma\left(\frac{3}{4}\right)}{8 \Gamma\left(\frac{7}{4}\right)}$

sympy.integrate(sympy.sin(x**2),(x, -sympy.oo, sympy.oo))

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

sympy.oo

$\displaystyle \infty$

sympy.limit(sympy.sin(x)/x, x, 0)

$\displaystyle 1$

sympy.limit(sympy.sin(x)/x, x, sympy.oo)

$\displaystyle 0$

sympy.Limit(sympy.sin(x)/x, x, 0)

$\displaystyle \lim_{x \to 0^+}\left(\frac{\sin{\left(x \right)}}{x}\right)$

sympy.solve(x**2-2)

[-sqrt(2), sqrt(2)]

y = sympy.Function('y')

type(y)

sympy.core.function.UndefinedFunction

t = sympy.symbols('t')

sympy.dsolve( sympy.Eq(y(t).diff(t,t) - y(t), sympy.exp(t)), y(t))

$\displaystyle y{\left(t \right)} = C_{2} e^{- t} + \left(C_{1} + \frac{t}{2}\right) e^{t}$

sympy.integrate(sympy.cos(x**2),x)

$\displaystyle \frac{\sqrt{2} \sqrt{\pi} C\left(\frac{\sqrt{2} x}{\sqrt{\pi}}\right) \Gamma\left(\frac{1}{4}\right)}{8 \Gamma\left(\frac{5}{4}\right)}$

sympy.Integrate(sympy.cos(x**2),x)

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

AttributeError                            Traceback (most recent call last)

<ipython-input-36-f1dbaae609e2> in <module>
----> 1 sympy.Integrate(sympy.cos(x**2),x)


AttributeError: module 'sympy' has no attribute 'Integrate'

sympy.Integral(sympy.cos(x**2),x)

$\displaystyle \int { \cos^{2}(x) }\, dx$

sympy.integrate(sympy.cos(x)**2,x)

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

sympy.latex(sympy.Integral(sympy.cos(x**2),x))

'\\int \\cos{\\left(x^{2} \\right)}\\, dx'

2.7 Calculus

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

sympy.init_printing(use_unicode=True)

sympy.diff(sympy.exp(x**3))

$\displaystyle 3 x^{2} e^{x^{3}}$

sympy.diff(x**4, x, 3)

$\displaystyle 24 x$

expr=sympy.exp(x*y*z)
expr

$\displaystyle e^{x y z}$

sympy.diff(expr, x,y,2,z,4)

$\displaystyle x^{3} y^{2} \left(x^{3} y^{3} z^{3} + 14 x^{2} y^{2} z^{2} + 52 x y z + 48\right) e^{x y z}$

也可以用點語法

expr.diff(x,y,y,z,4)

$\displaystyle x^{3} y^{2} \left(x^{3} y^{3} z^{3} + 14 x^{2} y^{2} z^{2} + 52 x y z + 48\right) e^{x y z}$

deriv = sympy.Derivative(expr, x,y,y,z,4)

deriv

$\displaystyle \frac{\partial^{7}}{\partial z^{4}\partial y^{2}\partial x} e^{x y z}$

deriv.doit()

$\displaystyle x^{3} y^{2} \left(x^{3} y^{3} z^{3} + 14 x^{2} y^{2} z^{2} + 52 x y z + 48\right) e^{x y z}$

畢文斌毛悅悅的書 Python 漫遊數學王國

這本書有較詳細的講解使用 python, SymPy, SciPy, SymPy, matplotlib 等
進行微積分, 線性代數, 統計機率, 以及 OR(運籌學, 作業研究)的各種例子, 有詳細的程式碼, 很有參考價值:

畢文斌毛悅悅編著, Python 漫遊數學王國, 清華大學出版,
繁體字版為: Python AI人員必修的科學計算, 深智出版, 202302.

Reference

SciPy 的入門說明講義不太好找, 在 https://scipy-lectures.org/index.html 這裡 link.
Scipy Lecture Notes, 1.6.5 Optimization and fit: scipy.optimize, https://scipy-lectures.org/index.html# link.
1.6.5.3. Finding the roots of a scalar function, https://scipy-lectures.org/intro/scipy.html#finding-the-roots-of-a-scalar-function link.
小狐狸事務所: Python 學習筆記 : SciPy 測試 (一) : 科學常數, http://yhhuang1966.blogspot.com/2020/04/python-scipy.html link.
以下此站主要介紹使用 sympy 算微積分的課程講義:
陳擎文教學網：python求解數學式（高中數學，大學數學，工程數學，微積分）
方程式，多項式，函數，三角函數，微分，積分，泰勒級數，傅立葉級數，傅立葉積分，傅立葉轉換，一階微分方程式，二階微分方程式，三階微分方程式，偏微分方程式，聯立方程式矩陣解法，張量，向量的坐標轉換，複數，極坐標
https://acupun.site/lecture/python_math/index.htm#chp5 link.
SymPy 的官網線上手冊: Welcome to SymPy’s documentation!, https://docs.sympy.org/latest/index.html link .
SymPy 的模組可以在此查看: SymPy Modules Reference, https://docs.sympy.org/latest/modules/index.html link
SymPy 延遲求值的方式,
Prevent expression evaluation, Ref: Advanced Expression Manipulation, https://docs.sympy.org/latest/tutorial/manipulation.html#prevent-expression-evaluation link
賴鵬仁, GeoGebra 與動態幾何 20200925 申請優良教材 link
賴鵬仁, 數學文書處理 1 link
賴鵬仁, 用免費的電腦資源協助數學的教學,學習與探索_復華中學教師營_中山大學應數系高中數學人才班_2021, https://blog.csdn.net/m0_47985483/article/details/113790840 link
GeoGebra 官網: https://www.geogebra.org/ link
intro-en_4_2_簡體中譯前半段_2015_103-GeoGebra-使用手冊之出處_有tangram, Geogebra 官網之 intro-en_4_2.pdf 之前半段, 由志工翻譯成中文, link
GeoGebra 指令完整列表 DivisorsList等等基隆女中上課教材1030319: link
geogebra与matlab,浅谈Geogebra在大学数学教学中的应用 link
畢文斌毛悅悅編著, Python 漫遊數學王國, 清華大學出版,
繁體字版為: Python AI人員必修的科學計算, 深智出版, 202302.

賴博伯

关注

0
点赞
踩
2

收藏

觉得还不错? 一键收藏
0
评论
用 Python+SciPy+SymPy+GeoGebra 執行微積分計算 1

GeoGebra 在中學數理科到大一微積分相當夠Python 中執行微積分或符號運算的主力是 SymPyPython 中執行數值計算的主力是 SciPy(含NumPy)如果考量到1. 免費開源社群穩定 2. 兼具數值與符號運算 3.又能延伸到研究所,建議同時使用 GeoGebra 搭配 Python+Scipy+SymPy, (或數值與符號運算無法兼具: R, Octave, Maxima 等, 如果不在意付費, 則 Matlab, Maple, Mathematica 搭配, 都可.)
复制链接

扫一扫