如何使用python計(jì)算常微分方程？

常用形式

創(chuàng)新互聯(lián)公司主營(yíng)平輿網(wǎng)站建設(shè)的網(wǎng)絡(luò)公司,主營(yíng)網(wǎng)站建設(shè)方案,重慶APP開發(fā),平輿h5小程序開發(fā)搭建,平輿網(wǎng)站營(yíng)銷推廣歡迎平輿等地區(qū)企業(yè)咨詢

odeint(func, y0, t,args,Dfun)

一般這種形式就夠用了。

下面是官方的例子，求解的是

D(D(y1))-t*y1=0

為了方便，采取D=d/dt。如果我們令初值

y1(0) = 1.0/3**(2.0/3.0)/gamma(2.0/3.0)

D(y1)(0) = -1.0/3**(1.0/3.0)/gamma(1.0/3.0)

這個(gè)微分方程的解y1=airy(t)。

令D(y1)=y0，就有這個(gè)常微分方程組。

D(y0)=t*y1

D(y1)=y0

Python求解該微分方程。

from scipy.integrate import odeint

from scipy.special import gamma, airy

y1_0 = 1.0/3**(2.0/3.0)/gamma(2.0/3.0)

y0_0 = -1.0/3**(1.0/3.0)/gamma(1.0/3.0)

y0 = [y0_0, y1_0]

def func(y, t):

... return [t*y[1],y[0]]

def gradient(y,t):

... return [[0,t],[1,0]]

x = arange(0,4.0, 0.01)

t = x

ychk = airy(x)[0]

y = odeint(func, y0, t)

y2 = odeint(func, y0, t, Dfun=gradient)

print ychk[:36:6]

[ 0.355028 0.339511 0.324068 0.308763 0.293658 0.278806]

print y[:36:6,1]

[ 0.355028 0.339511 0.324067 0.308763 0.293658 0.278806]

print y2[:36:6,1]

[ 0.355028 0.339511 0.324067 0.308763 0.293658 0.278806]

得到的解與精確值相比，誤差相當(dāng)小。

=======================================================================================================

args是額外的參數(shù)。

用法請(qǐng)參看下面的例子。這是一個(gè)洛侖茲曲線的求解，并且用matplotlib繪出空間曲線圖。（來(lái)自《python科學(xué)計(jì)算》）

from scipy.integrate import odeint

import numpy as np

def lorenz(w, t, p, r, b):

# 給出位置矢量w，和三個(gè)參數(shù)p, r, b 計(jì)算出

# dx/dt, dy/dt, dz/dt 的值

x, y, z = w

# 直接與lorenz 的計(jì)算公式對(duì)應(yīng)

return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])

t = np.arange(0, 30, 0.01) # 創(chuàng)建時(shí)間點(diǎn)

# 調(diào)用ode 對(duì)lorenz 進(jìn)行求解, 用兩個(gè)不同的初始值

track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))

track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0))

# 繪圖

from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt

fig = plt.figure()

ax = Axes3D(fig)

ax.plot(track1[:,0], track1[:,1], track1[:,2])

ax.plot(track2[:,0], track2[:,1], track2[:,2])

plt.show()

===========================================================================

scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)

計(jì)算常微分方程（組）

使用 FORTRAN庫(kù)odepack中的lsoda解常微分方程。這個(gè)函數(shù)一般求解初值問(wèn)題。

參數(shù)：

func : callable(y, t0, ...) 計(jì)算y在t0 處的導(dǎo)數(shù)。

y0 : 數(shù)組 y的初值條件（可以是矢量）

t : 數(shù)組 為求出y，這是一個(gè)時(shí)間點(diǎn)的序列。初值點(diǎn)應(yīng)該是這個(gè)序列的第一個(gè)元素。

args : 元組 func的額外參數(shù)

Dfun : callable(y, t0, ...) 函數(shù)的梯度(Jacobian)。即雅可比多項(xiàng)式。

col_deriv : boolean. True，Dfun定義列向?qū)?shù)（更快），否則Dfun會(huì)定義橫排導(dǎo)數(shù)

full_output : boolean 可選輸出，如果為True 則返回一個(gè)字典，作為第二輸出。

printmessg : boolean 是否打印convergence 消息。

返回： y : array, shape (len(y0), len(t))

數(shù)組，包含y值，每一個(gè)對(duì)應(yīng)于時(shí)間序列中的t。初值y0 在第一排。

infodict : 字典，只有full_output == True 時(shí)，才會(huì)返回。

字典包含額為的輸出信息。

鍵值：

‘hu’ vector of step sizes successfully used for each time step.

‘tcur’ vector with the value of t reached for each time step. (will always be at least as large as the input times).

‘tolsf’ vector of tolerance scale factors, greater than 1.0, computed when a request for too much accuracy was detected.

‘tsw’ value of t at the time of the last method switch (given for each time step)

‘nst’ cumulative number of time steps

‘nfe’ cumulative number of function evaluations for each time step

‘nje’ cumulative number of jacobian evaluations for each time step

‘nqu’ a vector of method orders for each successful step.

‘imxer’index of the component of largest magnitude in the weighted local error vector (e / ewt) on an error return, -1 otherwise.

‘lenrw’ the length of the double work array required.

‘leniw’ the length of integer work array required.

‘mused’a vector of method indicators for each successful time step: 1: adams (nonstiff), 2: bdf (stiff)

其他參數(shù)，官方網(wǎng)站和文檔都沒(méi)有明確說(shuō)明。相關(guān)的資料，暫時(shí)也找不到。

Python怎么做最優(yōu)化

一、概觀

scipy中的optimize子包中提供了常用的最優(yōu)化算法函數(shù)實(shí)現(xiàn)。我們可以直接調(diào)用這些函數(shù)完成我們的優(yōu)化問(wèn)題。optimize中函數(shù)最典型的特點(diǎn)就是能夠從函數(shù)名稱上看出是使用了什么算法。下面optimize包中函數(shù)的概覽：

1.非線性最優(yōu)化

fmin -- 簡(jiǎn)單Nelder-Mead算法

fmin_powell -- 改進(jìn)型Powell法

fmin_bfgs -- 擬Newton法

fmin_cg -- 非線性共軛梯度法

fmin_ncg -- 線性搜索Newton共軛梯度法

leastsq -- 最小二乘

2.有約束的多元函數(shù)問(wèn)題

fmin_l_bfgs_b ---使用L-BFGS-B算法

fmin_tnc ---梯度信息

fmin_cobyla ---線性逼近

fmin_slsqp ---序列最小二乘法

nnls ---解|| Ax - b ||_2 for x=0

3.全局優(yōu)化

anneal ---模擬退火算法

brute --強(qiáng)力法

4.標(biāo)量函數(shù)

fminbound

brent

golden

bracket

5.擬合

curve_fit-- 使用非線性最小二乘法擬合

6.標(biāo)量函數(shù)求根

brentq ---classic Brent (1973)

brenth ---A variation on the classic Brent（1980）ridder ---Ridder是提出這個(gè)算法的人名

bisect ---二分法

newton ---牛頓法

fixed_point

7.多維函數(shù)求根

fsolve ---通用

broyden1 ---Broyden’s first Jacobian approximation.

broyden2 ---Broyden’s second Jacobian approximationnewton_krylov ---Krylov approximation for inverse Jacobiananderson ---extended Anderson mixing

excitingmixing ---tuned diagonal Jacobian approximationlinearmixing ---scalar Jacobian approximationdiagbroyden ---diagonal Broyden Jacobian approximation8.實(shí)用函數(shù)

line_search ---找到滿足強(qiáng)Wolfe的alpha值

check_grad ---通過(guò)和前向有限差分逼近比較檢查梯度函數(shù)的正確性二、實(shí)戰(zhàn)非線性最優(yōu)化

fmin完整的調(diào)用形式是：

fmin(func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, full_output=0, disp=1, retall=0, callback=None)不過(guò)我們最常使用的就是前兩個(gè)參數(shù)。一個(gè)描述優(yōu)化問(wèn)題的函數(shù)以及初值。后面的那些參數(shù)我們也很容易理解。如果您能用到，請(qǐng)自己研究。下面研究一個(gè)最簡(jiǎn)單的問(wèn)題，來(lái)感受這個(gè)函數(shù)的使用方法：f(x)=x**2-4*x+8，我們知道，這個(gè)函數(shù)的最小值是4，在x=2的時(shí)候取到。

from scipy.optimize import fmin #引入優(yōu)化包def myfunc(x):

return x**2-4*x+8 #定義函數(shù)

x0 = [1.3] #猜一個(gè)初值

xopt = fmin(myfunc, x0) #求解

print xopt #打印結(jié)果

運(yùn)行之后，給出的結(jié)果是：

Optimization terminated successfully.

Current function value: 4.000000

Iterations: 16

Function evaluations: 32

[ 2.00001953]

程序準(zhǔn)確的計(jì)算得出了最小值，不過(guò)最小值點(diǎn)并不是嚴(yán)格的2，這應(yīng)該是由二進(jìn)制機(jī)器編碼誤差造成的。

除了fmin_ncg必須提供梯度信息外，其他幾個(gè)函數(shù)的調(diào)用大同小異，完全類似。我們不妨做一個(gè)對(duì)比：

from scipy.optimize import fmin,fmin_powell,fmin_bfgs,fmin_cgdef myfunc(x):

return x**2-4*x+8

x0 = [1.3]

xopt1 = fmin(myfunc, x0)

print xopt1

print

xopt2 = fmin_powell(myfunc, x0)

print xopt2

print

xopt3 = fmin_bfgs(myfunc, x0)

print xopt3

print

xopt4 = fmin_cg(myfunc,x0)

print xopt4

給出的結(jié)果是：

Optimization terminated successfully.

Current function value: 4.000000

Iterations: 16

Function evaluations: 32

[ 2.00001953]

Optimization terminated successfully.

Current function value: 4.000000

Iterations: 2

Function evaluations: 53

1.99999999997

Optimization terminated successfully.

Current function value: 4.000000

Iterations: 2

Function evaluations: 12

Gradient evaluations: 4

[ 2.00000001]

Optimization terminated successfully.

Current function value: 4.000000

Iterations: 2

Function evaluations: 15

Gradient evaluations: 5

[ 2.]

我們可以根據(jù)給出的消息直觀的判斷算法的執(zhí)行情況。每一種算法數(shù)學(xué)上的問(wèn)題，請(qǐng)自己看書學(xué)習(xí)。個(gè)人感覺，如果不是純研究數(shù)學(xué)的工作，沒(méi)必要搞清楚那些推導(dǎo)以及定理云云。不過(guò)，必須了解每一種算法的優(yōu)劣以及能力所及。在使用的時(shí)候，不妨多種算法都使用一下，看看效果分別如何，同時(shí)，還可以互相印證算法失效的問(wèn)題。

在from scipy.optimize import fmin之后，就可以使用help(fmin)來(lái)查看fmin的幫助信息了。幫助信息中沒(méi)有例子，但是給出了每一個(gè)參數(shù)的含義說(shuō)明，這是調(diào)用函數(shù)時(shí)候的最有價(jià)值參考。

有源碼研究癖好的，或者當(dāng)你需要改進(jìn)這些已經(jīng)實(shí)現(xiàn)的算法的時(shí)候，可能需要查看optimize中的每種算法的源代碼。在這里：https:/ / github. com/scipy/scipy/blob/master/scipy/optimize/optimize.py聰明的你肯定發(fā)現(xiàn)了，順著這個(gè)鏈接往上一級(jí)、再往上一級(jí)，你會(huì)找到scipy的幾乎所有源碼！

python如何實(shí)現(xiàn)求函數(shù)的在一個(gè)連續(xù)區(qū)間的最值？

如果函數(shù)是確定的，可以用導(dǎo)數(shù)的方法進(jìn)行計(jì)算，但是如果函數(shù)是不確定的，就需要用優(yōu)化的方法來(lái)處理了，比如常用的梯度上升法，模擬退火等，希望可以幫到你。