橢圓微分方程python

❶ python里怎麼樣求解微分方程

有很多大學生問我，學習python有什麼用呢？我說：你至少可以用來解微分方程，如下面的例子,就是解決微分方程：
y"+a*y'+b*y=0

代碼如下：

[python]view plain

• #y"+a*y'+b*y=0

• fromscipy.integrateimportodeint

• frompylabimport*

• defderiv(y,t):#返回值是y和y的導數組成的數組

• a=-2.0

• b=-0.1

• returnarray([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計算橢圓面積

答: 首先要知道橢圓的面積公式， S=pi*a*b。其中pi約等於3.14ab，分別表示橢圓的長短軸和短半軸。 Python代碼如下所示。

總的來說，這道題並不算很難，題主是對Python的相關語法不太了解嘛。

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

常用形式
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是額外的參數。
用法請參看下面的例子。這是一個洛侖茲曲線的求解，並且用matplotlib繪出空間曲線圖。（來自《python科學計算》）
from scipy.integrate import odeint
import numpy as np
def lorenz(w, t, p, r, b):
# 給出位置矢量w，和三個參數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) # 創建時間點
# 調用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解常微分方程。這個函數一般求解初值問題。

參數：

func : callable(y, t0, ...) 計算y在t0 處的導數。
y0 : 數組 y的初值條件（可以是矢量）
t : 數組 為求出y，這是一個時間點的序列。初值點應該是這個序列的第一個元素。
args : 元組 func的額外參數
Dfun : callable(y, t0, ...) 函數的梯度(Jacobian)。即雅可比多項式。
col_deriv : boolean. True，Dfun定義列向導數（更快），否則Dfun會定義橫排導數
full_output : boolean 可選輸出，如果為True 則返回一個字典，作為第二輸出。
printmessg : boolean 是否列印convergence 消息。

返回： y : array, shape (len(y0), len(t))
數組，包含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)
其他參數，官方網站和文檔都沒有明確說明。相關的資料，暫時也找不到。

❹ 常微分方程的解析解(方法歸納)以及基於Python的微分方程數值解算例實現

本文歸納常見常微分方程的解析解解法以及基於Python的微分方程數值解算例實現。

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

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

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

某些方程看似不可分離變數，但是經過換元之後，其實還是可分離變數的，不要被這種方程迷惑。

形如

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

解法： (直接套公式)

伯努利方程
形如

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


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

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

形如

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

求解此類方程分兩步：

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

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

解的情況分為以下三種：

情況一：方程有兩個不同的實數解

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

情況二：方程有一個二重解
假設該解等於 ，此時方程的通解是

情況三：方程有一對共軛復解
假設這對解是 , 此時方程的通解是

對於 和特徵根的情況，對特解的情況做如下歸納：

形如

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

求解此類方程分兩步：

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

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

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

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

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

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

求解得到兩個不同的實數特徵根: .

此時方程的齊次通解是

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

將上式代入控制方程有

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

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

於是,原方程的全解為

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

且有

將以上各式代入邊界條件

解此方程組可得: .

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

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

假定其解存在唯一,
為求解的近似值, 類似於前面的做法,

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

問題二:有限差分方法算出其數值解及誤差
對於 原問題 ,取步長 h=0.2 ,用 有限差分 求其 近似解 ,並將結果與 精確解y(x)=-x-1 進行比較.

因為

先以將區間劃分為5份為例，求出數值解

結果:

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

結果:

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

結果:

將區間劃分為 份, 即 時.

結果:

繪圖比較 時數值解與真解的差距：

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

結果：