建站
咨詢
售后
建站咨詢
新聞中心
python里怎么樣求解微分方程
有很多大學(xué)生問(wèn)我，學(xué)習(xí)python有什么用呢？我說(shuō)：你至少可以用來(lái)解微分方程，如下面的例子,就是解決微分方程：
站在用戶的角度思考問(wèn)題，與客戶深入溝通，找到勐臘網(wǎng)站設(shè)計(jì)與勐臘網(wǎng)站推廣的解決方案，憑借多年的經(jīng)驗(yàn)，讓設(shè)計(jì)與互聯(lián)網(wǎng)技術(shù)結(jié)合，創(chuàng)造個(gè)性化、用戶體驗(yàn)好的作品，建站類型包括：成都網(wǎng)站制作、成都網(wǎng)站設(shè)計(jì)、企業(yè)官網(wǎng)、英文網(wǎng)站、手機(jī)端網(wǎng)站、網(wǎng)站推廣、域名注冊(cè)、虛擬空間、企業(yè)郵箱。業(yè)務(wù)覆蓋勐臘地區(qū)。
y"+a*y'+b*y=0?
代碼如下：
[python]?view plain?copy
#y"+a*y'+b*y=0
from?scipy.integrate?import?odeint
from?pylab?import?*
def?deriv(y,t):????????#?返回值是y和y的導(dǎo)數(shù)組成的數(shù)組
a?=?-2.0
b?=?-0.1
return?array([?y[1],?a*y[0]+b*y[1]?])
time?=?linspace(0.0,50.0,1000)
yinit?=?array([0.0005,0.2])?????#?初值
y?=?odeint(deriv,yinit,time)
figure()
plot(time,y[:,0],label='y')????#y[:,0]即返回值的第一列，是y的值。label是為了顯示legend用的。
plot(time,y[:,1],label="y'")?????#y[:,1]即返回值的第二列，是y’的值
xlabel('t')
ylabel('y')
legend()
show()
輸出結(jié)果如下：
常微分方程的解析解(方法歸納)以及基于Python的微分方程數(shù)值解算例實(shí)現(xiàn)
本文歸納常見(jiàn)常微分方程的解析解解法以及基于Python的微分方程數(shù)值解算例實(shí)現(xiàn)。
考慮常微分方程的解析解法，我們一般可以將其歸納為如下幾類：
這類微分方程可以變形成如下形式：
兩邊同時(shí)積分即可解出函數(shù)，難點(diǎn)主要在于不定積分，是最簡(jiǎn)單的微分方程。
某些方程看似不可分離變量，但是經(jīng)過(guò)換元之后，其實(shí)還是可分離變量的，不要被這種方程迷惑。
形如
的方程叫做一階線性微分方程，若 為0，則方程齊次，否則稱為非齊次。
解法： (直接套公式)
伯努利方程
形如
的方程稱為伯努利方程，這種方程可以通過(guò)以下步驟化為一階線性微分方程：
令 , 方程兩邊同時(shí)乘以 ，得到
即 .
這就將伯努利方程歸結(jié)為可以套公式的一階線性微分方程。
形如
的方程稱為二階常系數(shù)微分方程，若 ，則方程稱為齊次的，反之稱為非齊次的。以下默認(rèn)方程是非齊次的。
求解此類方程分兩步：
原方程的解 = 齊次通解 + 非齊次特解
首先假設(shè) .用特征方程法，寫出對(duì)應(yīng)的特征方程并且求解：
解的情況分為以下三種：
情況一：方程有兩個(gè)不同的實(shí)數(shù)解
假設(shè)兩個(gè)實(shí)數(shù)解分別是 , 此時(shí)方程的通解是
情況二：方程有一個(gè)二重解
假設(shè)該解等于 ，此時(shí)方程的通解是
情況三：方程有一對(duì)共軛復(fù)解
假設(shè)這對(duì)解是 , 此時(shí)方程的通解是
對(duì)于 和特征根的情況，對(duì)特解的情況做如下歸納：
形如
的方程叫做高階常系數(shù)微分方程，若 ，則方程是齊次的，否則是非齊次的。下面默認(rèn)方程是非齊次的。
求解此類方程分兩步：
原方程的解 = 齊次通解 + 非齊次特解
考慮帶有第三類邊界條件的二階常系數(shù)微分方程邊值問(wèn)題
問(wèn)題一:兩點(diǎn)邊值問(wèn)題的解析解
由于此方程是非齊次的,故 求解此類方程分兩步：
原方程的解 = 齊次通解 + 非齊次特解
首先假設(shè) . 用特征方程法，寫出對(duì)應(yīng)的特征方程
求解得到兩個(gè)不同的實(shí)數(shù)特征根: .
此時(shí)方程的齊次通解是
由于 . 所以非齊次特解形式為
將上式代入控制方程有
求解得: , 即非齊次特解為 .
原方程的解 = 齊次通解 + 非齊次特解
于是,原方程的全解為
因?yàn)樵搯?wèn)題給出的是第三類邊界條件,故需要求解的導(dǎo)函數(shù)
且有
將以上各式代入邊界條件
解此方程組可得: .
綜上所述,原兩點(diǎn)邊值問(wèn)題的解為
對(duì)一般的二階微分方程邊值問(wèn)題
假定其解存在唯一,
為求解的近似值, 類似于前面的做法,
考慮帶有第三類邊界條件的二階常系數(shù)微分方程邊值問(wèn)題
問(wèn)題二:有限差分方法算出其數(shù)值解及誤差
對(duì)于 原問(wèn)題 ,取步長(zhǎng) h=0.2 ,用 有限差分 求其 近似解 ,并將結(jié)果與 精確解y(x)=-x-1 進(jìn)行比較.
因?yàn)?/p>
先以將區(qū)間劃分為5份為例，求出數(shù)值解
結(jié)果:
是不是解出數(shù)值解就完事了呢？當(dāng)然不是。我們可以將問(wèn)題的差分格式解與問(wèn)題的真解進(jìn)行比較，以得到解的可靠性。通過(guò)數(shù)學(xué)計(jì)算我們得到問(wèn)題的真解為 ，現(xiàn)用范數(shù)來(lái)衡量誤差的大小：
結(jié)果:
接下來(lái)繪圖比較 時(shí)數(shù)值解與真解的差距：
結(jié)果:
將區(qū)間劃分為 份, 即 時(shí).
結(jié)果:
繪圖比較 時(shí)數(shù)值解與真解的差距：
最后，我們還可以從數(shù)學(xué)的角度分析所采用的差分格式的一些性質(zhì)。因?yàn)椴罘指袷降恼`差為 ， 所以理論上來(lái)說(shuō)網(wǎng)格每加密一倍，與真解的誤差大致會(huì)縮小到原來(lái)的 . 下討論網(wǎng)格加密時(shí)的變化：
結(jié)果：
如何使用python計(jì)算常微分方程？
常用形式
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í)也找不到。
本文標(biāo)題：python函數(shù)微分 python微積分函數(shù)
分享鏈接：http://ef60e0e.cn/article/doipjjp.html
