python解常微分方程_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