摘要：在學(xué)習(xí)與科研中，經(jīng)常會遇到一些數(shù)學(xué)運算問題，使用計算機完成運算具有速度快和準(zhǔn)確性高的優(yōu)勢。Python的Numpy包具有強大的科學(xué)運算功能，且具有其他許多主流科學(xué)計算語言不具備的免費、開源、輕量級和靈活的特點。本文使用Python語言的NumPy庫，解決數(shù)學(xué)運算問題中的線性方程組問題、積分問題、微分問題及矩陣化簡問題，結(jié)果準(zhǔn)確快捷，具有一定的借鑒意義。

1.Sympy庫簡介

SymPy一個用于符號型數(shù)學(xué)計算（symbolic mathematics）的Python庫。它旨在成為一個功能齊全的計算機代數(shù)系統(tǒng)（Computer Algebra System，CAS），同時保持代碼簡潔、易于理解和擴展。SymPy完全是用Python寫的，并不需要外部的庫。

本文選擇Sympy庫的原因在于：

• 免費：該庫基于BSD開源許可，免費且開源；
• 基于Python：該庫完全是用Python寫就，并以Python作為該庫操作語言；
• 輕量級：為了使Sympy簡單易用，該庫僅基于mpmath庫（一個純Python庫，用于浮點運算）；
• 靈活性：除了用作交互工具，還可插入其他應(yīng)用或軟件拓展功能中。

具體說來，如果x與y未曾賦值，那么下列語句就會報錯

#測試語句
print(x+y)

而符號運算則不同，符號運算多用于公式推導(dǎo)，不需要賦值，此時使用Sympy進行符號運算便具有方便快捷的優(yōu)勢，如下述語句便不會報錯。

#測試語句
x=Symbol('x')
y=Symbol('y')
print(x+y)

2 SymPy庫解決數(shù)學(xué)運算問題實現(xiàn)

2.1 求解線性方程組

解方程的功能主要是使用Sympy中solve函數(shù)實現(xiàn)。以式(1)為例，求解過程如下：

（1） 符號表示

SymPy庫中使用Symbol函數(shù)定義符號變量，

from sympy import *
x=Symbol('x')
y=Symbol('y')
 
#或者用如下語句
x,y=Symbol('x y')#第二個用空格隔開

（2）方程表示

使用代碼表示數(shù)學(xué)符號與手寫體的數(shù)學(xué)運算符號存在一定的差異，下面列舉常用的運算符：

• 加號
• 加號 +
• 減號 -
• 除號 /
• 乘號 *
• 指數(shù) **
• 對數(shù) log()
• e的指數(shù)次冪 exp()

對于長的表達(dá)式，如果不確定運算符的優(yōu)先級，可以加入小括號提升其優(yōu)先級。由于需要將表達(dá)式都轉(zhuǎn)化成右端等于0,這里把常數(shù)3和7移到等式左邊。題目中表達(dá)式可表示為：

2*x-y-3=0
3*x+y-7=0

(3)使用Solve函數(shù)解方程

在使用Solve函數(shù)解方程之前，我們先來看一下Solve函數(shù)的定義。Solve函數(shù)的第一個參數(shù)是要解的方程，要求右端等于0，第二個參數(shù)是未知數(shù)。

對于式(1)的求解，代碼如下：

from sympy import *
x = Symbol('x')
y = Symbol('y')
print(solve([2*x-y-3,3*x+y-7],[x,y]))

2.2 求解微積分問題

2.2.1 求解極限問題

在2.1中通過一個簡單的二元一次方程組求解熟悉了該庫求解數(shù)學(xué)問題的基本過程，下面本文通過示例，講解使用SymPy庫求解微積分的過程。

求解式（2）所示的極限問題，需要用到limit函數(shù)求極限。

(1)符號及方程表示

引入Sympy庫并定義n為符號變量與2.1中一致。

from sympy import *
n = Symbol('n')
s = ((n+3)/(n+2))**n

(2)利用limit函數(shù)求極限

首先我們介紹limit函數(shù)的調(diào)用格式：limit(e, z, z0, dir='+'),e為任意表達(dá)式，表示求取e(z)在點z0處的極限，dir='+'表示取右極限，die='-'則表示取左極限。則上式的求解代碼可表示如下：

from sympy import *
n = Symbol('n')
s = ((n+3)/(n+2))**n

print(limit(s,n,oo)) #無窮的表示方法是兩個小寫的字母o

2.2.2 求解定積分

(1)符號表示

from sympy import *
t = Symbol(t)
x = Symbol(x)

(2)方程表示

m = integrate(sin(t)/(pi-t),(t,0,x))
n = integrate(m,(x,0,pi)) 

完整代碼如下：

from sympy import *
t = Symbol(t)
x = Symbol(x)
m = integrate(sin(t)/(pi-t),(t,0,x))
n = integrate(m,(x,0,pi))
print(n)

2.2.3 求解微分問題

如求取的通解

（1）符號表示

這里與之前不同的是增加了函數(shù)的表示(用f(x)表示y)，即例題中的y還有微分表示

from sympy import *
f = Function('f')
x = Symbol('x')

y'的表示方法由以下代碼組成

diff(f(x),x)

這里對diff函數(shù)稍作介紹：

上面是求一階導(dǎo)的方法，求解高階導(dǎo)的方法如下所示：

>>> diff(x**3,x)
3*x**2
>>> diff(x**3,x,1)
3*x**2
>>> diff(x**3,x,2)
6*x
>>> diff(x**3,x,3)
6
>>> diff(x**3,x,4)
0

即改變第三個參數(shù)即可。

下面繼續(xù)我們的解題過程。

#左端
diff(f(x),x)
#看一下
print(diff(f(x),x))
 
#result
#d
#--(f(x))
#dx
 
#右端
2*f(x)*x

用dsolve函數(shù)解微分方程

dsolve函數(shù)是用來解決微分方程（differential equation）的函數(shù)。

函數(shù)的一個用法為：

dsolve(eq, f(x))

第一個參數(shù)為微分方程（要先將等式移項為右端為0的形式)。第二個參數(shù)為要解的函數(shù)(在微分方程中)

舉個例子:

>>> from sympy import *
>>> f = Function('f')
>>> x = Symbol('x')
>>> pprint(2*x-diff(f(x),x))
   d
2*x - --(f(x))
   dx
>>> dsolve(2*x - diff(f(x),x), f(x))
#result
#Eq(f(x), C1 + x**2)

這樣，我們可以將我們要解的題目，用以下代碼表示。

dsolve(diff(f(x),x) - 2*f(x)*x, f(x))

結(jié)果為：

Eq(f(x), C1*exp(x**2))
#即f(x) = C1*exp(x**2)

對比答案可以發(fā)現(xiàn)正確。

完整代碼:

from sympy import *
f = function('f')
x = Symbol('x')
print(dsolve(diff(f(x),x)-2*f(x)*x,f(x))

2.2.4 矩陣化簡

平時線性代數(shù)問題中我們會遇到化簡問題，雖然不算難，但著實麻煩。而且，出一點錯就會導(dǎo)致

結(jié)果出錯。不過好運的是SymPy提供了相關(guān)的支持。

例題：

符號表示與矩陣表示

from sympy import *
x1,x2,x3 = symbols('x1 x2 x3')
a11,a12,a13,a22,a23,a33 = symbols('a11 a12 a13 a22 a23 a33')
m = Matrix([[x1,x2,x3]])
n = Matrix([[a11,a12,a13],[a12,a22,a23],[a13,a23,a33]])
v = Matrix([[x1],[x2],[x3]])

注意m的表示，需要有兩個中括號

化簡實現(xiàn)

print(m*n*v)

得到的是：

Matrix([[x1*(a11*x1 + a12*x2 + a13*x3) + x2*(a12*x1 + a22*x2 + a23*x3) + x3*(a13*x1 + a23*x2 + a33*x3)]])

使用

f = m * n * v
print f[0]

可以進一步得到化簡后的式子

也許你要問我要化簡后在計算怎么辦？下面我就舉個例子。

如果上式中x1，x2，x3均等于1，則可這樣代入。

from sympy import *
x1,x2,x3 = symbols('x1 x2 x3')
a11,a12,a13,a22,a23,a33 = symbols('a11 a12 a13 a22 a23 a33')
m = Matrix([[x1, x2, x3]])
n = Matrix([[a11, a12, a13], [a12, a22, a23], [a13, a23, a33]])
v = Matrix([[x1], [x2], [x3]])
f = m * n * v
print f[0].subs({x1:1, x2:1, x3:1})

可得

a11 + 2*a12 + 2*a13 + a22 + 2*a23 + a33

以上就是本文的全部內(nèi)容，希望對大家的學(xué)習(xí)有所幫助，也希望大家多多支持腳本之家。

最新評論