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

1.Sympy庫簡(jiǎn)介

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

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

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

具體說來，如果x與y未曾賦值，那么下列語句就會(huì)報(bào)錯(cuò)

#測(cè)試語句
print(x+y)

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

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

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

2.1 求解線性方程組

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

（1） 符號(hào)表示

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

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

（2）方程表示

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

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

對(duì)于長(zhǎng)的表達(dá)式，如果不確定運(yùn)算符的優(yōu)先級(jí)，可以加入小括號(hào)提升其優(yōu)先級(jí)。由于需要將表達(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ù)的第一個(gè)參數(shù)是要解的方程，要求右端等于0，第二個(gè)參數(shù)是未知數(shù)。

對(duì)于式(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中通過一個(gè)簡(jiǎn)單的二元一次方程組求解熟悉了該庫求解數(shù)學(xué)問題的基本過程，下面本文通過示例，講解使用SymPy庫求解微積分的過程。

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

(1)符號(hào)及方程表示

引入Sympy庫并定義n為符號(hào)變量與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)在點(diǎn)z0處的極限，dir='+'表示取右極限，die='-'則表示取左極限。則上式的求解代碼可表示如下：

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

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

2.2.2 求解定積分

(1)符號(hào)表示

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）符號(hào)表示

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

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

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

diff(f(x),x)

這里對(duì)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

即改變第三個(gè)參數(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ù)的一個(gè)用法為：

dsolve(eq, f(x))

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

舉個(gè)例子:

>>> 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)

對(duì)比答案可以發(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 矩陣化簡(jiǎn)

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

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

例題：

符號(hào)表示與矩陣表示

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的表示，需要有兩個(gè)中括號(hào)

化簡(jiǎn)實(shí)現(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]

可以進(jìn)一步得到化簡(jiǎn)后的式子

也許你要問我要化簡(jiǎn)后在計(jì)算怎么辦？下面我就舉個(gè)例子。

如果上式中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

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

最新評(píng)論