Python3 SciPy解常微分方程 用Matplotlib演示

Python科學計算 簡單記錄幾篇筆記 用SciPy解常微分方程，在jupyter notebook用matplotlib演示，如下須要注意的幾點：

• integrate模塊提供的odeint函數
• Anaconda 3的jupyter notebook上
• matplotlib 2D 繪製求解 牛頓冷卻定律
• matplotlib 3D 繪製求解 洛倫茲吸引子

jupyter notebook上繪製2D圖表

import numpy as np
import matplotlib.pyplot as plt

x1 = np.linspace(0.0, 5.0)
x2 = np.linspace(0.0, 2.0)

y1 = np.cos(2 * np.pi * x1) * np.exp(-x1)
y2 = np.cos(2 * np.pi * x2)

plt.subplot(2, 1, 1)
plt.plot(x1, y1, '.-')
plt.title('A tale of 2 subplots')
plt.ylabel('Damped oscillation')

plt.subplot(2, 1, 2)
plt.plot(x2, y2, '-')
plt.xlabel('time (s)')
plt.ylabel('Undamped')

plt.show()

牛頓冷卻定律做爲最基礎的微分方程，就用它來演示第一個例子

T′(t)=−α(T(t)−H)T(t)=H+(T0−H)e−αtT′(t)=−α(T(t)−H)T(t)=H+(T0−H)e−αt

# -*- coding:utf-8 -*-

from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np

from IPython import display

# 冷卻定律的微分方程
def cooling_law_equ(w, t, a, H):
    return -1 * a * (w - H)

# 冷卻定律求解獲得的溫度temp關於時間t的函數
def cooling_law_func(t, a, H, T0):
    return H + (T0 - H) * np.e ** (-a * t)

t = np.arange(0, 10, 0.01)
initial_temp = (90) #初始溫度
temp = odeint(cooling_law_equ, initial_temp, t, args=(0.5, 30)) #冷卻係數和環境溫度
temp1 = cooling_law_func(t, 0.5, 30, initial_temp) #推導的函數與scipy計算的結果對比
plt.subplot(2, 1, 1)
plt.plot(t, temp)
plt.ylabel("temperature")

plt.subplot(2, 1, 2)
plt.plot(t, temp1)
plt.xlabel("time")
plt.ylabel("temperature")

display.Latex("牛頓冷卻定律 $T'(t)=-a(T(t)- H)$)（上）和 $T(t)=H+(T_0-H)e^{-at}$（下）")
plt.show()

 

jupyter notebook上繪製3D圖表

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt

mpl.rcParams['legend.fontsize'] = 10

fig = plt.figure()
ax = fig.gca(projection='3d')
theta = np.linspace(-4 * np.pi, 4 * np.pi, 100)
z = np.linspace(-2, 2, 100)
r = z**2 + 1
x = r * np.sin(theta)
y = r * np.cos(theta)
ax.plot(x, y, z, label='parametric curve')
ax.legend()

plt.show()

 

三維空間下洛倫茲吸引子的演示

dxdt=σ⋅(y−x)dydt=x⋅(ρ−z)−ydzdt=xy−βzdxdt=σ⋅(y−x)dydt=x⋅(ρ−z)−ydzdt=xy−βz 

%matplotlib inline 
from scipy.integrate import odeint
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt

from IPython import display

fig = plt.figure()
ax = fig.gca(projection='3d')

def lorenz(w, t, p, r, b):
    # 位置矢量w， 三個參數p, r, b
    x, y ,z = w.tolist()
    # 分別計算dx/dt, dy/dt, dz/dt
    return p * (y-x), x*(r-z)-y, x*y-b*z

t = np.arange(0, 30, 0.02)
initial_val = (0.0, 1.00, 0.0)

track = odeint(lorenz, initial_val, t, args=(10.0, 28.0, 3.0))
X, Y, Z = track[:,0], track[:,1], track[:,2]
ax.plot(X, Y, Z, label='lorenz')
ax.legend()

display.Latex(r"$\frac{dx}{dt}=\sigma\cdot(y-x) \\ \frac{dy}{dt}=x\cdot(\rho-z)-y \\ \frac{dz}{dt}=xy-\beta z$")

 

版權聲明：本文爲博主原創文章，未經博主容許不得轉載。 https://blog.csdn.net/m454078356/article/details/77776092

文章標籤： python

我的分類： Python

相關熱詞： python3∑ python3 python3與 python3包 python3庫