python應用-scipy,numpy,sympy計算微積分

python應用-scipy,numpy,sympy計算微積分

今天來講一下使用python進行微積分運算，python有很多科學計算庫都可以進行微積分運算，當然如果知曉微積分計算的原理也可以自己編程實現。

下面我們用三種方式進行積分運算圓周率pi

numpy計算pi

import os
import numpy as np
#pi=4(1-1/3+1/5-1/7+1/9-.......)

n = 100000 
print(np.sum(4.0 / np.r_[1:n:4, -3:-n:-4]))
#3.141572693

講解一下上面的代碼
首先，這里的編程思路來源於一個公司，代碼中也有注解

pi=4((1-1/3+1/5-1/7+1/9-1/11+1/13…)
這個公式怎么來的，實話說我也不知道
上面的代碼思路就是，不可能說按照上面的公式計算無數項從1-1/3+1/5…到無窮項，所以我們就盡可能地去計算，去逼近真實結果，事實上，別的微積分計算函數也不能說就能計算無窮多項，其實也只是一種盡可能地逼近。上面的代碼是取了1-1/3…1/99999。取到這個范圍。
還要提到一個方法np.r_[]是指按行連接兩個矩陣

所以如果我們呢自己去編程其實也是這樣的一個思路

下面我給大家舉一個很典型的例子。

計算微函數 $y=x^2$從0-5的積分，積分思路：將0-5的 $y=x^2$圖像進行切分，切成1000000個小矩形再求他們的面積和，說到這里很多人可能想起了高等數學教材上有關於這一節的積分知識。

下面給出代碼：

import os
import numpy as np


#pi=4(1-1/3+1/5-1/7+1/9-.......)

n=1000000
print(sum((5/1000000)*(np.linspace(0,5,1000000)**2)))

如果你用數學積分的知識計算，積分結果應該為125/3，上述結果的精確度已經達到0.00001了，下面給出圖形並給出繪圖代碼方便大家形象化了解

x=np.linspace(0,5,1000)
y=x**2
plt.figure()
plt.xlabel('x')
plt.ylabel('y')
plt.plot(x,y,label="$y=x^2$")

plt.legend()
x=np.linspace(0,5,100)

for i in x:
    x_z=[]
    y_z=[]
    x_z.append(i)
    y_z.append(0) 
    x_z.append(i)
    y_z.append(i**2)
    plt.plot(x_z,y_z)

plt.show()

運行結果：

下面介紹的積分方法都是直接調用方法了

在下面的例子中，使用SciPy中提供的數值積分函數quad()計算π：

$\pi=\int_{-1}^1\sqrt{1-x^2}dx$

scipy積分代碼：

from scipy.integrate import quad 
print(quad(lambda x:(1-x**2)**0.5, -1, 1)[0] * 2 )
#3.141592653589797

如上面代碼，quad，接受的第一個參數為一個函數，或者說一個積分公式，后面兩個參數分別為積分下限和上限

下面用SymPy提供的符號積分函數integrate()對上面的公式進行積分運算，可以看到 運算的結果為符號表示的π：
代碼：

from sympy import symbols, integrate, sqrt 
x = symbols("x") 
print(integrate(sqrt(1-x**2), (x, -1, 1)) * 2)

積分的方法為intergrate，我們需要先定義一個變量，然后傳遞參數，第一個參數為積分公式，第二個參數為積分上下限。