如何使用python計算常微分方程？

常用形式

公司主營業(yè)務：網(wǎng)站設計、做網(wǎng)站、移動網(wǎng)站開發(fā)等業(yè)務。幫助企業(yè)客戶真正實現(xiàn)互聯(lián)網(wǎng)宣傳，提高企業(yè)的競爭能力。創(chuàng)新互聯(lián)是一支青春激揚、勤奮敬業(yè)、活力青春激揚、勤奮敬業(yè)、活力澎湃、和諧高效的團隊。公司秉承以“開放、自由、嚴謹、自律”為核心的企業(yè)文化，感謝他們對我們的高要求，感謝他們從不同領域給我們帶來的挑戰(zhàn)，讓我們激情的團隊有機會用頭腦與智慧不斷的給客戶帶來驚喜。創(chuàng)新互聯(lián)推出山南免費做網(wǎng)站回饋大家。

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)

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

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

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]

得到的解與精確值相比，誤差相當小。

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

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

用法請參看下面的例子。這是一個洛侖茲曲線的求解，并且用matplotlib繪出空間曲線圖。（來自《python科學計算》）

from scipy.integrate import odeint

import numpy as np

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

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

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

x, y, z = w

# 直接與lorenz 的計算公式對應

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

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

# 調用ode 對lorenz 進行求解, 用兩個不同的初始值

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)

計算常微分方程（組）

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

參數(shù)：

func : callable(y, t0, ...) 計算y在t0 處的導數(shù)。

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

t : 數(shù)組 為求出y，這是一個時間點的序列。初值點應該是這個序列的第一個元素。

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

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

col_deriv : boolean. True，Dfun定義列向導數(shù)（更快），否則Dfun會定義橫排導數(shù)

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

printmessg : boolean 是否打印convergence 消息。

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

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

infodict : 字典，只有full_output == True 時，才會返回。

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

鍵值：

‘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)站和文檔都沒有明確說明。相關的資料，暫時也找不到。

python里怎么樣求解微分方程

有很多大學生問我，學習python有什么用呢？我說：你至少可以用來解微分方程，如下面的例子,就是解決微分方程：

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的導數(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()

輸出結果如下：

常微分方程的解析解(方法歸納)以及基于Python的微分方程數(shù)值解算例實現(xiàn)

本文歸納常見常微分方程的解析解解法以及基于Python的微分方程數(shù)值解算例實現(xiàn)。

考慮常微分方程的解析解法，我們一般可以將其歸納為如下幾類：

這類微分方程可以變形成如下形式：

兩邊同時積分即可解出函數(shù)，難點主要在于不定積分，是最簡單的微分方程。

某些方程看似不可分離變量，但是經(jīng)過換元之后，其實還是可分離變量的，不要被這種方程迷惑。

形如

的方程叫做一階線性微分方程，若 為0，則方程齊次，否則稱為非齊次。

解法： (直接套公式)

伯努利方程

形如

的方程稱為伯努利方程，這種方程可以通過以下步驟化為一階線性微分方程：

令 , 方程兩邊同時乘以 ，得到

即 .

這就將伯努利方程歸結為可以套公式的一階線性微分方程。

形如

的方程稱為二階常系數(shù)微分方程，若 ，則方程稱為齊次的，反之稱為非齊次的。以下默認方程是非齊次的。

求解此類方程分兩步：

原方程的解 = 齊次通解 + 非齊次特解

首先假設 .用特征方程法，寫出對應的特征方程并且求解：

解的情況分為以下三種：

情況一：方程有兩個不同的實數(shù)解

假設兩個實數(shù)解分別是 , 此時方程的通解是

情況二：方程有一個二重解

假設該解等于 ，此時方程的通解是

情況三：方程有一對共軛復解

假設這對解是 , 此時方程的通解是

對于 和特征根的情況，對特解的情況做如下歸納：

形如

的方程叫做高階常系數(shù)微分方程，若 ，則方程是齊次的，否則是非齊次的。下面默認方程是非齊次的。

求解此類方程分兩步：

原方程的解 = 齊次通解 + 非齊次特解

考慮帶有第三類邊界條件的二階常系數(shù)微分方程邊值問題

問題一:兩點邊值問題的解析解

由于此方程是非齊次的,故 求解此類方程分兩步：

原方程的解 = 齊次通解 + 非齊次特解

首先假設 . 用特征方程法，寫出對應的特征方程

求解得到兩個不同的實數(shù)特征根: .

此時方程的齊次通解是

由于 . 所以非齊次特解形式為

將上式代入控制方程有

求解得: , 即非齊次特解為 .

原方程的解 = 齊次通解 + 非齊次特解

于是,原方程的全解為

因為該問題給出的是第三類邊界條件,故需要求解的導函數(shù)

且有

將以上各式代入邊界條件

解此方程組可得: .

綜上所述,原兩點邊值問題的解為

對一般的二階微分方程邊值問題

假定其解存在唯一,

為求解的近似值, 類似于前面的做法,

考慮帶有第三類邊界條件的二階常系數(shù)微分方程邊值問題

問題二:有限差分方法算出其數(shù)值解及誤差

對于 原問題 ,取步長 h=0.2 ,用 有限差分 求其 近似解 ,并將結果與 精確解y(x)=-x-1 進行比較.

因為

先以將區(qū)間劃分為5份為例，求出數(shù)值解

結果:

是不是解出數(shù)值解就完事了呢？當然不是。我們可以將問題的差分格式解與問題的真解進行比較，以得到解的可靠性。通過數(shù)學計算我們得到問題的真解為 ，現(xiàn)用范數(shù)來衡量誤差的大小：

結果:

接下來繪圖比較 時數(shù)值解與真解的差距：

結果:

將區(qū)間劃分為 份, 即 時.

結果:

繪圖比較 時數(shù)值解與真解的差距：

最后，我們還可以從數(shù)學的角度分析所采用的差分格式的一些性質。因為差分格式的誤差為 ， 所以理論上來說網(wǎng)格每加密一倍，與真解的誤差大致會縮小到原來的 . 下討論網(wǎng)格加密時的變化：

結果：

如何用python求導數(shù)

打開python運行環(huán)境。

導入微分的模塊包：from sympy import *。

定義符號變量：x = symbols('x')

定義一個函數(shù)：f = x**9

diff = diff(f,x)求導

最后輸入diff，即可顯示其變量值了。

眾多python培訓視頻，盡在python學習網(wǎng)，歡迎在線學習！

python 的scipy 里的 odeint 這個求微分方程的函數(shù)怎么用啊

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)

實際使用中，還是主要使用前三個參數(shù)，即微分方程的描寫函數(shù)、初值和需要求解函數(shù)值對應的的時間點。接收數(shù)組形式。這個函數(shù)，要求微分方程必須化為標準形式，即dy/dt=f(y,t,)。

fromscipyimportodeint

y=odeint(dy/dt=r*y*(1-y/k),y(0)=0.1,t)

對于微分方程全還給老師了，

網(wǎng)站名稱：python的求函數(shù)微分,python計算微分方程
當前網(wǎng)址：http://chinadenli.net/article6/dsijgig.html

成都網(wǎng)站建設公司_創(chuàng)新互聯(lián)，為您提供品牌網(wǎng)站設計、外貿建站、做網(wǎng)站、企業(yè)建站、電子商務、建站公司

廣告

聲明：本網(wǎng)站發(fā)布的內容（圖片、視頻和文字）以用戶投稿、用戶轉載內容為主，如果涉及侵權請盡快告知，我們將會在第一時間刪除。文章觀點不代表本網(wǎng)站立場，如需處理請聯(lián)系客服。電話：028-86922220；郵箱：631063699@qq.com。內容未經(jīng)允許不得轉載，或轉載時需注明來源： 創(chuàng)新互聯(lián)