Python在數(shù)值分析中的應(yīng)用

摘 要：數(shù)值分析是理科專(zhuān)業(yè)重要的基礎(chǔ)課之一，Python語(yǔ)言在數(shù)值分析中的應(yīng)用廣泛而深入。NumPy是Python中最基礎(chǔ)的科學(xué)計(jì)算庫(kù)之一，它提供了多維數(shù)組對(duì)象，具備許多數(shù)學(xué)函數(shù)，能夠處理大型矩陣。SciPy是建立在NumPy之上的庫(kù)，在工程計(jì)算方面有豐富的功能子模塊。Matplotlib是一款用于數(shù)據(jù)可視化的Python繪圖庫(kù)，它支持多種圖表，包括線圖、散點(diǎn)圖、柱狀圖、餅圖、等高線圖等。本文給出了應(yīng)用Python程序解決數(shù)值分析常見(jiàn)問(wèn)題的示例，展示了這幾個(gè)流行庫(kù)的基本功能，為數(shù)值分析知識(shí)的學(xué)習(xí)和掌握提供了借鑒。

關(guān)鍵詞：Python語(yǔ)言；數(shù)值分析；數(shù)值解；解析解

數(shù)值分析和Python都是人工智能和數(shù)據(jù)科學(xué)等許多理工科專(zhuān)業(yè)的必修科目。數(shù)值分析也稱(chēng)為計(jì)算數(shù)學(xué)，是研究如何在計(jì)算機(jī)上求解問(wèn)題的一個(gè)數(shù)學(xué)分支。Python是目前最流行的語(yǔ)言之一，它簡(jiǎn)單易學(xué)、免費(fèi)開(kāi)源、各類(lèi)模塊資源非常豐富。隨著Python的興起，早期使用其他語(yǔ)言（如Fortran和MATLAB）做數(shù)值計(jì)算的許多用戶(hù)開(kāi)始使用Python語(yǔ)言。Python在數(shù)值分析方面的模塊有很多，功能十分強(qiáng)大，用戶(hù)不必關(guān)注當(dāng)中的函數(shù)是如何實(shí)現(xiàn)的，可直接調(diào)用。以下通過(guò)幾個(gè)實(shí)例，介紹Python程序在數(shù)值分析中的應(yīng)用。

1 解線性方程組

線性方程組有很多種解法，如高斯消去法和雅可比迭代法［1］。使用NumPy庫(kù)可以直接求解，而不必考慮具體的求解過(guò)程。函數(shù)numpy.linalg.solve（）用于解決形如Ax=B的線性方程組，其中A為系數(shù)矩陣，B為等號(hào)右側(cè)數(shù)值組成的列向量，x為未知數(shù)列向量。例如，有如下方程組：

5x1-4x2+x3=2

-4x1+6x2-4x3=-1

x1-4x2+6x3=-1

相應(yīng)的解該方程組的Python代碼如下：

import numpy as np

A=np.array（［［5，-4，1］，［-4，6，-4］，［1，-4，6］］）

B=np.array（［2，-1，-1］）

x=np.linalg.solve（A，B）

print（x）

運(yùn)行結(jié)果為：

［0.33333333-0.16666667-0.33333333］

即x1=1/3，x2=-1/6，x3=-1/3，代入原方程驗(yàn)證，這個(gè)結(jié)果是正確的。

2 數(shù)值微分

有的函數(shù)導(dǎo)數(shù)沒(méi)有解析解，則函數(shù)的差商可以作為其導(dǎo)數(shù)的近似值，即：

y′=df（x）dx≈f（x+h）-f（x）h

式中的h必須是一個(gè)很小的數(shù)，當(dāng)h無(wú)限趨近于0時(shí)，上式就是導(dǎo)數(shù)的定義式。上式是向前差商，也可以采用向后差商或者中心差商近似代替導(dǎo)數(shù)［2］。

下面的代碼中，自定義了差商函數(shù)derivative（），參數(shù)f和a表示求函數(shù)f在a處的差商，參數(shù)md的默認(rèn)值是central，表示中心差商，md也可以取forward（向前差商）或者backward（向后差商）。我們知道正弦函數(shù)的導(dǎo)數(shù)是余弦函數(shù)，即，sin'（x）=cos（x），下面的代碼使用差商函數(shù)求sin'（x）的數(shù)值解：

def derivative（f，a，md='central'，h=0.01）：

if md=='central'：

return（f（a+h）-f（a-h））/（2*h）

elif md=='forward'：

return（f（a+h）-f（a））/h

else：#'backward'：

return（f（a）-f（a-h））/h

x=np.linspace（0，5*np.pi，100）

dy_dx=derivative（np.sin，x）

y=np.cos（x）

plt.plot（x，y，\"b.\"，label='y=cos（x）'）

plt.plot（x，dy_dx，label=\"y=sin'（x）\"）

plt.legend（）

plt.grid（True）

plt.show（）

代碼的后半部分使用子模塊matplotlib.pyplot這個(gè)最常用的Python繪圖工具將cos（x）曲線和sin'（x）數(shù)值解的散點(diǎn)畫(huà)在了同一幅圖中，以便于對(duì)比。如圖1所示，連續(xù)的曲線表示cos（x），散點(diǎn)*表示sin'（x）數(shù)值解，顯然，sin'（x）數(shù)值解與cos（x）非常吻合。

3 數(shù)值積分

積分運(yùn)算，包括牛頓柯特斯公式、復(fù)合（梯形與辛普森）求積公式和龍貝格公式等［3］。例如，求解下面的定積分值：

∫2π0xsin（x）dx

SciPy庫(kù)的子模塊integrate主要用于數(shù)值積分的運(yùn)算。將函數(shù)表達(dá)式f（x），積分下限a與積分上限b代入函數(shù)integrate.quad（f，a，b）當(dāng)中，可計(jì)算出定積分值。這里，f（x）=xsin（x），a=0，b=2π，相應(yīng)的Python代碼如下：

import math

from scipy import integrate

def f（x）：

return x*math.sin（x）

Ifx，err=integrate.quad（f，0，2*math.pi）

print（'定積分為：'，Ifx）

print（'誤差為：'，err）

運(yùn)行結(jié)果為：

定積分為：-6.2831853071795845

誤差為：1.382265978143112e-13

看來(lái)誤差非常小，說(shuō)明計(jì)算精度比較高。同時(shí)不難看出，定積分的值約為-6.28，十分接近-2π，實(shí)際上，通過(guò)分步積分容易求得，其準(zhǔn)確值就是-2π。SymPy是一個(gè)以符號(hào)運(yùn)算為主的Python庫(kù)，支持解析解運(yùn)算，應(yīng)用它對(duì)前面的定積分進(jìn)行符號(hào)運(yùn)算求解，相應(yīng)的代碼如下：

from sympy import *

x=Symbol（\"x\"）

Fx=x*sin（x）

IFx =integrate（Fx，（x，0，2*pi））

print（IFx）

運(yùn)行結(jié)果為：

-2*pi

這里的pi是SymPy中的符號(hào)，表示π。

4 解常微分方程

下面是一個(gè)有解析解的、并有初值條件的一階微分方程：y′=x2+x-y，y（0）=0。

將方程變形為：y′+y-x2-x=0。

相應(yīng)求解析解的Python程序如下，代碼中的f1就是上式的等號(hào)左邊的表達(dá)式，ics為初值條件y（0）=0：

from sympy import *

x=symbols（\"x\"）

y=symbols（\"y\"，cls=Function）

f1=diff（y（x），x）+y（x）-x**2-x

result=dsolve（f1，y（x），ics={y（0）：0}）

print（result）

運(yùn)行結(jié)果為：

Eq（y（x），x**2-x+1-exp（-x））

即解析解為：

y=x2-x+1-e-x

多數(shù)微分方程是沒(méi)有解析解的，只能求數(shù)值解，方法有歐拉法、梯形法和龍格庫(kù)塔法等［4］。Python子模塊scipy.integrate中的函數(shù)odeint（）用于計(jì)算微分方程的數(shù)值解。

上述原始微分方程的等號(hào)右邊為x2+x-y，也就是y′的表達(dá)式，將其定義為函數(shù)dfun（y，x），x的散點(diǎn)取0到1.1，步長(zhǎng)為0.05，相應(yīng)的求該微分方程數(shù)值解的代碼如下：

from scipy.integrate import odeint

def dfun（y，x）：

return x*x+x-y

x=np.arange（0，1.1，0.05）

f2=x*x-x+1-np.exp（-x）

ans=odeint（dfun，0，x）

plt.plot（x，ans，\"b*\"，label='Numerical'）

plt.plot（x，f2，\"red\"，label='Analytical'）

plt.grid（）

plt.legend（loc=\"best\"）

plt.show（）

為了比較數(shù)值解和解析解的相近程度，在上面的代碼里包含了解析解（用f2表示），使用子模塊matplotlib.pyplot在同一張圖中畫(huà)出了數(shù)值解的散點(diǎn)和解析解的曲線，如圖2所示。

可見(jiàn)，函數(shù)odeint（）的計(jì)算結(jié)果是比較準(zhǔn)確的。

5 最小二乘擬合

最小二乘擬合的基本原理是，通過(guò)給定數(shù)據(jù)點(diǎn)，尋找一個(gè)函數(shù)的近似表達(dá)式。具體地說(shuō)，有一組（m個(gè)）實(shí)驗(yàn)數(shù)據(jù)（xi，yi），假設(shè)它們之間應(yīng)該滿(mǎn)足某函數(shù)關(guān)系y=f（x），例如，如果是線性關(guān)系，則f（x）=kx+b，k和b就是待定參數(shù)。找到一組k和b的值，使得下面的誤差的平方和最?。?/p>

s（k，b）=∑mi=1［yi-f（xi）］2

有了k和b的值就得到了函數(shù)的近似表達(dá)式，這種算法就稱(chēng)為最小二乘擬合［5］。

下面的程序，先使用隨機(jī)函數(shù)，生成101個(gè)近似地處在一條直線的散點(diǎn)，然后利用scipy.optimize模塊中的函數(shù)curve_fit（）找到最佳參數(shù)值：

x=np.linspace（0，1，101）

y=1+x+0.1*np.random.random（len（x））

def func（x，k，b）：

return k*x+b

k，b=curve_fit（func，xdata=x，ydata=y）［0］

print（'k='，k）

print（'b='，b）

plt.plot（x，y，\"b.\"）

plt.plot（x，k*x+b，\"r\"）

plt.grid（）

plt.show（）

運(yùn)行結(jié)果為：

k=1.0013726989903353

b=1.0448271910479061

在同一張圖中畫(huà)出了模擬的散點(diǎn)圖和擬合出來(lái)的直線，如圖3所示?？梢?jiàn)，散點(diǎn)在直線兩側(cè)的分布比較均勻，說(shuō)明擬合效果較好。

6 特征值與特征向量

求矩陣特征值和特征向量的基本方法有：冪法與反冪法、QR算法以及雅可比法（針對(duì)實(shí)對(duì)稱(chēng)矩陣）。子模塊numpy.linalg中的函數(shù)eig（）用于求解特征值和特征向量。例如，有如下矩陣A：

A=120

2-11

013

相應(yīng)的求特征值和特征向量的Python程序如下：

import numpy as np

from numpy.linalg import eig

A=np.array（［［1，2，0］，［2，-1，1］，［0，1，3］］）

eigenvalue，eigenvector=eig（A）

print（eigenvalue）

print（eigenvector）

運(yùn)行結(jié)果為：

［-2.37228132 2. 3.37228132］

［［0.50369186 0.81649658 0.28218405］

［-0.84929533 0.40824829 0.33470998］

［0.1580884 -0.40824829 0.89907808］］

可以驗(yàn)證程序運(yùn)行結(jié)果是正確的。

結(jié)語(yǔ)

以上例子都比較簡(jiǎn)單，實(shí)際上，使用Python還可以進(jìn)行插值運(yùn)算、傅里葉變換、求解偏微分方程和非線性方程組等，在數(shù)值積分方面，還可以計(jì)算多重積分。本文未提及的數(shù)值計(jì)算方面的Python（子）模塊還有很多，作為數(shù)值分析的學(xué)習(xí)者和工作者，不能僅關(guān)心如何使用Python相應(yīng)的模塊工具求解，更應(yīng)該掌握每種解法背后的數(shù)學(xué)原理及求解過(guò)程的來(lái)龍去脈，這樣才能確保在學(xué)習(xí)的道路上步伐邁得穩(wěn)且走得遠(yuǎn)。

參考文獻(xiàn)：

［1］李慶揚(yáng).數(shù)值分析［M］.5版.北京：清華大學(xué)出版社，2008：142188.

［2］RAJAT SHARMA.Numerical Methods in Python［EB/OL］.［20240704］.https：//medium.com/.

［3］朱曉臨.數(shù)值分析［M］.2版.合肥：中國(guó)科學(xué)技術(shù)大學(xué)出版社，2014：230248.

［4］孔慶凱.Python編程與數(shù)值方法［M］.北京：機(jī)械工業(yè)出版社，2023：289298.

［5］尹紅麗.Python科學(xué)計(jì)算數(shù)據(jù)處理與分析［M］.北京：人民郵電出版社，2023：237240.

作者簡(jiǎn)介：石也牧（2004— ），女，漢族，北京人，本科，研究方向：人工智能。