Python在數值分析中的應用

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

關鍵詞：Python語言；數值分析；數值解；解析解

數值分析和Python都是人工智能和數據科學等許多理工科專業的必修科目。數值分析也稱為計算數學，是研究如何在計算機上求解問題的一個數學分支。Python是目前最流行的語言之一，它簡單易學、免費開源、各類模塊資源非常豐富。隨著Python的興起，早期使用其他語言（如Fortran和MATLAB）做數值計算的許多用戶開始使用Python語言。Python在數值分析方面的模塊有很多，功能十分強大，用戶不必關注當中的函數是如何實現的，可直接調用。以下通過幾個實例，介紹Python程序在數值分析中的應用。

1 解線性方程組

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

5x1-4x2+x3=2

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

x1-4x2+6x3=-1

相應的解該方程組的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）

運行結果為：

［0.33333333-0.16666667-0.33333333］

即x1=1/3，x2=-1/6，x3=-1/3，代入原方程驗證，這個結果是正確的。

2 數值微分

有的函數導數沒有解析解，則函數的差商可以作為其導數的近似值，即：

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

式中的h必須是一個很小的數，當h無限趨近于0時，上式就是導數的定義式。上式是向前差商，也可以采用向后差商或者中心差商近似代替導數［2］。

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

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這個最常用的Python繪圖工具將cos（x）曲線和sin'（x）數值解的散點畫在了同一幅圖中，以便于對比。如圖1所示，連續的曲線表示cos（x），散點*表示sin'（x）數值解，顯然，sin'（x）數值解與cos（x）非常吻合。

3 數值積分

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

∫2π0xsin（x）dx

SciPy庫的子模塊integrate主要用于數值積分的運算。將函數表達式f（x），積分下限a與積分上限b代入函數integrate.quad（f，a，b）當中，可計算出定積分值。這里，f（x）=xsin（x），a=0，b=2π，相應的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）

運行結果為：

定積分為：-6.2831853071795845

誤差為：1.382265978143112e-13

看來誤差非常小，說明計算精度比較高。同時不難看出，定積分的值約為-6.28，十分接近-2π，實際上，通過分步積分容易求得，其準確值就是-2π。SymPy是一個以符號運算為主的Python庫，支持解析解運算，應用它對前面的定積分進行符號運算求解，相應的代碼如下：

from sympy import *

x=Symbol（\"x\"）

Fx=x*sin（x）

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

print（IFx）

運行結果為：

-2*pi

這里的pi是SymPy中的符號，表示π。

4 解常微分方程

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

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

相應求解析解的Python程序如下，代碼中的f1就是上式的等號左邊的表達式，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）

運行結果為：

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

即解析解為：

y=x2-x+1-e-x

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

上述原始微分方程的等號右邊為x2+x-y，也就是y′的表達式，將其定義為函數dfun（y，x），x的散點取0到1.1，步長為0.05，相應的求該微分方程數值解的代碼如下：

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（）

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

可見，函數odeint（）的計算結果是比較準確的。

5 最小二乘擬合

最小二乘擬合的基本原理是，通過給定數據點，尋找一個函數的近似表達式。具體地說，有一組（m個）實驗數據（xi，yi），假設它們之間應該滿足某函數關系y=f（x），例如，如果是線性關系，則f（x）=kx+b，k和b就是待定參數。找到一組k和b的值，使得下面的誤差的平方和最小：

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

有了k和b的值就得到了函數的近似表達式，這種算法就稱為最小二乘擬合［5］。

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

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（）

運行結果為：

k=1.0013726989903353

b=1.0448271910479061

在同一張圖中畫出了模擬的散點圖和擬合出來的直線，如圖3所示。可見，散點在直線兩側的分布比較均勻，說明擬合效果較好。

6 特征值與特征向量

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

A=120

2-11

013

相應的求特征值和特征向量的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）

運行結果為：

［-2.37228132 2. 3.37228132］

［［0.50369186 0.81649658 0.28218405］

［-0.84929533 0.40824829 0.33470998］

［0.1580884 -0.40824829 0.89907808］］

可以驗證程序運行結果是正確的。

結語

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

參考文獻：

［1］李慶揚.數值分析［M］.5版.北京：清華大學出版社，2008：142188.

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

［3］朱曉臨.數值分析［M］.2版.合肥：中國科學技術大學出版社，2014：230248.

［4］孔慶凱.Python編程與數值方法［M］.北京：機械工業出版社，2023：289298.

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

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