数据结构
网络
关系型数据库管理系统
操作系统
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 次查看
