資料結構
網路
關係型資料庫管理系統
作業系統
Java
MS Excel
iOS
HTML
CSS
Android
Python
C 語言程式設計
C++
C#
MongoDB
MySQL
Javascript
PHP
物理學
化學
生物學
數學
英語
經濟學
心理學
社會學
服裝設計
法律研究使用 Python 建模泰勒表方法
泰勒表方法是一種非常高效且優雅的方法，用於在考慮特定模板大小的情況下，獲得特定導數的有限差分格式。要理解它，必須非常清楚什麼是模板。
假設想要評估 $\mathrm{\frac{d^{2}f}{dx^{2}}}$，那麼在有限差分方法中，起點是泰勒級數。請參考下圖以更好地理解該方法。
點 $\mathrm{x_{i} \: + \: h}$ 處的泰勒級數展開式為
$$\mathrm{f(x_{i} \: + \: h) \: = \: f(x_{i}) \: + \: hf'(x_{i}) \: + \: \frac{h^{2}}{2!}f''(x_{i}) \: + \:\frac{h^{3}}{3!}f'''(x_{i}) \: + \:\frac{h^{3}}{3!}f'''(x_{i}) \: + \: \dotso}$$
$$\mathrm{f(x_{i} \: − \: h) \: = \: f(x_{i}) \: − \: hf'(x_{i}) \: + \: \frac{h^{2}}{2!}f''(x_{i}) \: − \:\frac{h^{3}}{3!}f'''(x_{i}) \: + \:\frac{h^{3}}{3!}f'''(x_{i}) \: − \: \dotso}$$
讓我們將以上兩個方程相加，並將包含二階導數的項分離出來
$$\mathrm{f''(x_{i}) \: = \: \frac{f(x_{i} \: + \: h) \: − \: 2f(x_{i}) \: + \: f(x_{i} \: − \: h)}{h^{2}}}$$
現在，這是考慮點 $\mathrm{f(x_{i} \: + \: h) \: , \: f(x_{i})}$ 和 $\mathrm{f(x_{i} \: − \: h)}$ 的二階導數近似值。這就是我們所說的模板。基本上，模板是用於數值逼近導數的一組點。為了簡化書寫任務，每當需要寫 $\mathrm{f(x_{i})}$ 時，都將簡單地寫成 $\mathrm{f_{i}}$，類似地，將 $\mathrm{f(x_{i} \: − \: h)}$ 寫成 $\mathrm{f_{i − 1}}$ 等等。
現在這很容易，但當任務是考慮模板中非常多的點數，並且必須針對此模板獲得某個導數的數值表示式時，問題就會變得非常複雜和困難。假設已經給出需要獲得三階導數 $\mathrm{(f'''(x_{i}))}$，並考慮模板 $\mathrm{f_{i} \: , \: f_{i+1} \: , \: f_{i+2} \: , \: f_{i+3}}$。那麼任務將變得非常繁瑣且耗時。因此，為了簡化任務，泰勒表方法變得非常方便。
該方法可以很容易地解釋如下
首先，將導數寫成如下形式
$$\mathrm{f'''(x_{i}) \: = \: a \: f_{i} \: + \: b \: f_{i+1} \: + \: c \: f_{i+2} \: + \: d \: f_{i+3}}$$
現在任務是找到 a、b、c 和 d。我們將導數及其係數排列成如下表格的形式。
用黃色突出顯示的單元格將填充對應於 $\mathrm{1^{st}}$ 列中顯示的模板的泰勒級數的係數。在紅色突出顯示的單元格中，我們將放置 0，除了對應於 $\mathrm{f'''(x_{i})}$ 的那個單元格。
第一行乘以 a，第二行乘以 b，第三行乘以 c，第四行乘以 d，然後全部相加並與最後一列（考慮維度）相等，得到以下方程
$$\mathrm{a \: + \: b \: + \: c \: + \: d \: = \: 0}$$
$$\mathrm{(0 \: + \: b \: + \: 2c \: + \: 2d) \: = \: 0/h}$$
$$\mathrm{(0 \: + \: b/2 \: + \: 2c \: + \: 9d/2) \: = \: 0/h^{\wedge}2}$$
$$\mathrm{(0 \: + \: b/6 \: + \: 4c/3 \: + \: 9d/2) \: = \: 0/h^{\wedge}3}$$
現在這些可以很容易地排列成矩陣形式，很容易求解。但仍然這是一個非常繁瑣的任務，如果模板大小進一步增加，那麼它將變得非常困難。因此，本文使用 Python 程式設計來簡化任務。
Python 中泰勒表的實現
我們將要做的第一件事是編寫一個程式碼來獲得泰勒級數的係數。通常，泰勒級數寫成。
$$\mathrm{f(x_{i} \: + \: \beta h) \: = \: f(x_{i}) \: + \: \beta f'(x_{i})h \: + \: \frac{\beta^{2}}{2!}f''(x_{i})h^{2} \: + \: \frac{\beta^{3}}{3!}f'''(x_{i})h^{3} \: + \: \frac{\beta^{4}}{4!}f''''(x_{i})h^{4} \: + \: \dots}$$
因此，上述級數的 $\mathrm{i^{th}}$ 係數可以寫成：$\mathrm{(\frac{\alpha^{i}}{i!})}$。所以首先我們將為它編寫函式，如下所示
Importing modules
from numpy import *
from math import *
# Taylor series function
def TSC(n,α):
"""
n: Number of terms in the series
α: coefficient of h
"""
# Empty array of coefficients
c=empty(n)
for i in range(n):
c[i]=α**i/factorial(i)
return c
函式 TSC 將接受級數中的項數和 $\mathrm{\alpha}$。
然後，我們將使用此函式來獲得不同模板元素的係數，即之前用黃色突出顯示的泰勒表的行。
為了演示它，假設想要使用模板 $\mathrm{(f_{i-1} \: , \: f_{i}, \: and \: f_{i+1})}$ 來評估 $\mathrm{d^{2}f/dx^{2}}$。我們將形成一個模板陣列，使其包含以下程式碼中提到的元素
$$\mathrm{f_{i} \: : \: 0}$$
$$\mathrm{f_{i − 1} \: : \: −1}$$
$$\mathrm{f_{i+1} \: : \: +1}$$
$$\mathrm{f_{i−2} \: : \: −2}$$
$$\mathrm{f_{i+3} \: : \: 3}$$
我們還將提供必須透過泰勒表逼近的導數的階數。然後，使用 for 迴圈填充表中的元素，如下所示
# Instruction for setting stencil #================================# # 1. f_i --> 0 # 2. f_(i-1)-->-1 # 3. f_(i+1)--> 1 # 4. f_(i-2)-->-2 # 5. f_(i+2)--> 2 #================================# # Forming Stencil array st=array([-1,0,1]) # Order of derivative to be evaluated order=2 # Number of rows and columns in the Taylor table nt=len(st) Taylor_Table=empty((nt,nt)) for i in range(nt): Taylor_Table[i,:]=TSC(nt,st[i]) Taylor_Table
輸出將如下所示
array([[1. , -1. , 0.5],
[1. , 0. , 0. ],
[1. , 1. , 0.5]])
現在必須如下設定向量 b
# Right hand vector b=zeros(nt) # Setting b as per derivatiev given b[order]=1 b
這將返回以下輸出
$$\mathrm{array([0. \: , \: 0. \: , \: 1.])}$$
在設定向量時，我們已確保給定導數正下方的那個值應為 1。但泰勒表矩陣和向量不能按原樣求解，因為它們沒有采用正確的形式，因此必須在應用任何內建矩陣求解方法之前進行轉置，然後如下求解
a=linalg.solve(transpose(Taylor_Table),transpose(b))
這將產生以下輸出
$$\mathrm{array([−0.5 \: , \: 0. \: , \: 0.5])}$$
但是，在最終答案中，我們應該得到這些值除以 h 的某些適當冪。這是藉助以下程式碼完成的
for i in range(nt):
print(f"a[{i}] = {a[i]}/h^{order}")
因此，最終答案將如下所示
a[0] = 1.0/h^2 a[1] = -2.0/h^2 a[2] = 1.0/h^2
因此，導數可以寫成
$$\mathrm{f''(x_{i}) \: = \: a[0]f_{i-1} \: + \: a[1]f_{i} \: + \: a[2]f_{i+1}}$$
$$\mathrm{f''(x_{i}) \: = \: \frac{f(x_{i} \: + \: h) \: − \: 2f(x_{i}) \: + \: f(x_{i} \: − \: h)}{h^{2}}}$$
完整程式碼如下
# Importing modules
from numpy import *
from math import *
# Taylor series function
#================================#
def TSC(n,α):
"""
n: Number of terms in the series
α: coefficient of h
"""
# Empty array of coefficients
c=empty(n)
for i in range(n):
c[i]=α**i/factorial(i)
return c
#================================#
# Instruction for setting stencil
# 1. f_i --> 0
# 2. f_(i-1)-->-1
# 3. f_(i+1)--> 1
# 4. f_(i-2)-->-2
# 5. F1_(i+2)--> 2
#================================#
# Forming the Stencil array
st=array([-1,0,1])
# Order of derivative term
order=2
# Number of rows and columns in the taylor table
nt=len(st)
Taylor_Table=empty((nt,nt))
for i in range(nt):
Taylor_Table[i,:]=TSC(nt,st[i])
Taylor_Table
#================================#
# Right hand vector
b=zeros(nt)
# Setting b as per derivatiev given
b[order]=1
#================================#
# Solving the matrix
a=linalg.solve(transpose(Taylor_Table),transpose(b))
#================================#
# Printing the final answer
for i in range(nt):
print(f"a[{i}] = {a[i]}/h^{order}")
#================================#
在上面的程式碼中，只需要更改第 33 行和第 36 行的模板和階數。
讓我們解決引言部分中提到的問題。
$$\mathrm{f'''(x_{i}) \: = \: a[0]f_{i} \: + \: a[1]f_{i+1} \: + \: a[2]f_{i+2} \: + \: a[3]f_{i+3}} $$
這裡，模板將為：st=array([0,1,2,3])，階數將為 3。程式輸出將如下所示
a[0] = -1.0/h^3 a[1] = 3.0/h^3 a[2] = -3.0/h^3 a[3] = 1.0/h^3
結論
在本教程中，使用 Python 程式設計實現了泰勒表方法。此外，首先解釋了該方法，然後解釋了 Python 策略。
208 次檢視
