当前位置： 首页 > news >正文

用Python进行数学建模（二）

news 来源：原创 2024/5/19 18:23:03

一、微分方程模型

微分方程是描述系统的状态随时间和空间演化的数学工具。物理中许多涉及变力的运动学、动力学问题，如空气的阻力为速度函数的落体运动等问题，很多可以用微分方程求解。微分方程在化学、工程学、经济学和人口统计等领域也有广泛应用。

具体来说，微分方程是指含有未知函数及其导数的关系式。

微分方程的数学建模其实并不复杂，基本过程就是分析题目属于哪一类问题、可以选择什么微分方程模型，然后如何使用现有的微分方程模型建模。

1.微分方程的数值解法

微分方程的数值求解是先把时间和空间离散化，然后将微分化为差分，建立递推关系，然后反复进行迭代计算，得到任意时间和空间的值。
在这里插入图片描述

2.SciPy 求解常微分方程（组）

1.一阶常微分方程（组）模型
在这里插入图片描述
2.scipy.integrate.odeint() 函数

函数原型
在这里插入图片描述

参数解释
func:callable(y, t, …) or callable(t, y, …)
计算y在t点的导数。如果签名是可调用的(t, y，…)，那么参数tfirst必须设置为True。
y0:array
y的初始条件 (可以是vector).
t:array
要解出y的时间点序列。初值点应该是这个序列的第一个元素。这个序列必须是单调递增或单调递减的;允许重复值。
args:tuple, optional
向导数函数 func 传递参数。当导数函数 f ( y , t , p 1 , p 2 , . . ) f(y,t,p1,p2,…)f(y,t,p1,p2,…) 包括可变参数 p1,p2… 时，通过 args =(p1,p2,…) 可以将参数p1,p2… 传递给导数函数 func。

odeint 的主要返回值
y: array
数组，形状为 (len(t),len(y0)，给出时间序列 t 中每个时刻的 y 值。

3.Scipy 求解一阶常微分方程

1.scipy.integrate.odeint() 求解常微分方程初值问题的步骤
在这里插入图片描述
2.求微分方程的数值解
在这里插入图片描述

from scipy.integrate import odeint  # 导入 scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as plt

def dy_dt(y, t):  # 定义函数 f(y,t)
    return np.sin(t**2)

y0 = [1]  # y0 = 1 也可以
t = np.arange(-10,10,0.01)  # (start,stop,step)
y = odeint(dy_dt, y0, t)  # 求解微分方程初值问题

# 绘图
plt.plot(t, y)
plt.title("scipy.integrate.odeint")
plt.show()

在这里插入图片描述

4.Scipy 求解一阶常微分方程组

1.求洛伦兹（Lorenz）方程的数值解
在这里插入图片描述
2.洛伦兹（Lorenz）方程问题的编程步骤
在这里插入图片描述
3.代码

from scipy.integrate import odeint  # 导入 scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


# 导数函数, 求 W=[x,y,z] 点的导数 dW/dt
def lorenz(W, t, p, r, b):
    x, y, z = W  # W=[x,y,z]
    dx_dt = p * (y - x)  # dx/dt = p*(y-x), p: sigma
    dy_dt = x * (r - z) - y  # dy/dt = x*(r-z)-y, r:rho
    dz_dt = x * y - b * z  # dz/dt = x*y - b*z, b;beta
    return np.array([dx_dt, dy_dt, dz_dt])


t = np.arange(0, 30, 0.01)  # 创建时间点 (start,stop,step)
paras = (10.0, 28.0, 3.0)  # 设置 Lorenz 方程中的参数 (p,r,b)

# 调用ode对lorenz进行求解, 用两个不同的初始值 W1、W2 分别求解
W1 = (0.0, 1.00, 0.0)  # 定义初值为 W1
track1 = odeint(lorenz, W1, t, args=(10.0, 28.0, 3.0))  # args 设置导数函数的参数
W2 = (0.0, 1.01, 0.0)  # 定义初值为 W2
track2 = odeint(lorenz, W2, t, args=paras)  # 通过 paras 传递导数函数的参数

# 绘图
fig = plt.figure()
ax = Axes3D(fig, auto_add_to_figure=False)
fig.add_axes(ax)
ax.plot(track1[:, 0], track1[:, 1], track1[:, 2], color='magenta')  # 绘制轨迹 1
ax.plot(track2[:, 0], track2[:, 1], track2[:, 2], color='deepskyblue')  # 绘制轨迹 2
ax.set_title("Lorenz attractor by scipy.integrate.odeint")
plt.show()

4.结果
在这里插入图片描述
在这里插入图片描述

5.Scipy 求解高阶常微分方程

高阶常微分方程，必须做变量替换，化为一阶微分方程组，再用 odeint 求数值解。
1.求二阶 RLC 振荡电路的数值解
在这里插入图片描述
2.二阶微分方程问题的编程步骤
在这里插入图片描述
3.代码

from scipy.integrate import odeint  # 导入 scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as plt


# 导数函数，求 Y=[u,v] 点的导数 dY/dt
def deriv(Y, t, a, w):
    u, v = Y  # Y=[u,v]
    dY_dt = [v, -2 * a * v - w * w * u]
    return dY_dt


t = np.arange(0, 20, 0.01)  # 创建时间点 (start,stop,step)
# 设置导数函数中的参数 (a, w)
paras1 = (1, 0.6)  # 过阻尼：a^2 - w^2 > 0
paras2 = (1, 1)  # 临界阻尼：a^2 - w^2 = 0
paras3 = (0.3, 1)  # 欠阻尼：a^2 - w^2 < 0

# 调用ode对进行求解, 用两个不同的初始值 W1、W2 分别求解
Y0 = (1.0, 0.0)  # 定义初值为 Y0=[u0,v0]
Y1 = odeint(deriv, Y0, t, args=paras1)  # args 设置导数函数的参数
Y2 = odeint(deriv, Y0, t, args=paras2)  # args 设置导数函数的参数
Y3 = odeint(deriv, Y0, t, args=paras3)  # args 设置导数函数的参数

# 绘图
plt.plot(t, Y1[:, 0], 'r-', label='u1(t)')
plt.plot(t, Y2[:, 0], 'b-', label='u2(t)')
plt.plot(t, Y3[:, 0], 'g-', label='u3(t)')
plt.plot(t, Y1[:, 1], 'r:', label='v1(t)')
plt.plot(t, Y2[:, 1], 'b:', label='v2(t)')
plt.plot(t, Y3[:, 1], 'g:', label='v3(t)')
plt.axis([0, 20, -0.8, 1.2])
plt.legend(loc='best')
plt.title("Second ODE by scipy.integrate.odeint")
plt.show()

4.结果
在这里插入图片描述
在这里插入图片描述

二、微分方程边值问题（BVP）

微分方程是指含有未知函数及其导数的关系式。

微分方程是描述系统的状态随时间和空间演化的数学工具。

微分方程分为初值问题和边值问题。初值问题是已知微分方程的初始条件，即自变量为零时的函数值，一般可以用欧拉法、龙哥库塔法来求解。边值问题则是已知微分方程的边界条件，即自变量在边界点时的函数值。

1.常微分方程边值问题的数学模型

在这里插入图片描述

2.SciPy 求解常微分方程边值问题

1.BVP 问题的标准形式
在这里插入图片描述
2.scipy.integrate.solve_bvp() 函数
可以求解一阶微分方程（组）的两点边值问题（第一类边界条件）。
在这里插入图片描述
solve_bvp 的主要参数：
func: callable fun(x, y, …) 　　导数函数 f ( y , x ) f(y,x)f(y,x) ， y 在 x 处的导数，以函数的形式表示。可以带有参数 p。
bc: callable bc(ya, yb, …) 　　边界条件，y 在两点边界的函数，以函数的形式表示。可以带有参数 p。
x: array：　　初始网格的序列，shape (m,)。必须是单调递增的实数序列，起止于两点边界值 xa，xb。
y: array：　　网格节点处函数值的初值，shape (n,m)，第 i 列对应于 x[i]。

solve_bvp 的主要返回值：
sol: PPoly 　　通过 PPoly （如三次连续样条函数）插值求出网格节点处的 y 值。
x: array 　　数组，形状为 (m,)，最终输出的网格节点。
y: array 　　二维数组，形状为 (n,m)，输出的网格节点处的 y 值。
yp: array 　　二维数组，形状为 (n,m)，输出的网格节点处的 y’ 值。

3.一阶常微分方程边值问题

1.标准化
在这里插入图片描述
2.编程步骤
在这里插入图片描述
3.代码

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


# 1. 求解微分方程组边值问题，DEMO
# y'' + abs(y) = 0, y(0)=0.5, y(4)=-1.5

# 导数函数，计算导数 dY/dx
def dydx(x, y):
    dy0 = y[1]
    dy1 = -abs(y[0])
    return np.vstack((dy0, dy1))


# 计算 边界条件
def boundCond(ya, yb):
    fa = 0.5  # 边界条件 y(xa=0) = 0.5
    fb = -1.5  # 边界条件 y(xb=4) = -1.5
    return np.array([ya[0] - fa, yb[0] - fb])


xa, xb = 0, 4  # 边界点 (xa,xb)
# fa, fb = 0.5, -1.5  # 边界点的 y值
xini = np.linspace(xa, xb, 11)  # 确定 x 的初值
yini = np.zeros((2, xini.size))  # 确定 y 的初值
res = solve_bvp(dydx, boundCond, xini, yini)  # 求解 BVP

xSol = np.linspace(xa, xb, 100)  # 输出的网格节点
ySol = res.sol(xSol)[0]  # 网格节点处的 y 值

plt.plot(xSol, ySol, label='y')
plt.xlabel("x")
plt.ylabel("y")
plt.title("scipy.integrate.solve_bvp")
plt.show()

4.执行结果
在这里插入图片描述

4.水滴横截面的形状

在这里插入图片描述

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


# 3. 求解微分方程边值问题，水滴的横截面
# 导数函数，计算 h=[h0,h1] 点的导数 dh/dx
def dhdx(x, h):
    # 计算 dh0/dx, dh1/dx 的值
    dh0 = h[1]  # 计算 dh0/dx
    dh1 = (h[0] - 1) * (1 + h[1] * h[1]) ** 1.5  # 计算 dh1/dx
    return np.vstack((dh0, dh1))


# 计算 边界条件
def boundCond(ha, hb):
    # ha = 0  # 边界条件：h0(x=-1) = 0
    # hb = 0  # 边界条件：h0(x=1) = 0
    return np.array([ha[0], hb[0]])


xa, xb = -1, 1  # 边界点 (xa=0, xb=1)
xini = np.linspace(xa, xb, 11)  # 设置 x 的初值
hini = np.zeros((2, xini.size))  # 设置 h 的初值

res = solve_bvp(dhdx, boundCond, xini, hini)  # 求解 BVP
# scipy.integrate.solve_bvp(fun, bc, x, y,..)
#   fun(x, y, ..), 导数函数 f(y,x)，y在 x 处的导数。
#   bc(ya, yb, ..), 边界条件，y 在两点边界的函数。
#   x: shape (m)，初始网格的序列，起止于两点边界值 xa，xb。
#   y: shape (n,m)，网格节点处函数值的初值，第 i 列对应于 x[i]。

xSol = np.linspace(xa, xb, 100)  # 输出的网格节点
hSol = res.sol(xSol)[0]  # 网格节点处的 h 值
plt.plot(xSol, hSol, label='h(x)')
plt.xlabel("x")
plt.ylabel("h(x)")
plt.axis([-1, 1, 0, 1])
plt.title("Cross section of water drop by BVP")
plt.show()

在这里插入图片描述

5.带有未知参数的微分方程边值问题

在这里插入图片描述
在这里插入图片描述

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


# 4. 求解微分方程组边值问题，Mathieu 方程
# y0' = y1, y1' = -(lam-2*q*cos(2x))y0)
# y0(0)=1, y1(0)=0, y1(pi)=0

# 导数函数，计算导数 dY/dx
def dydx(x, y, p):  # p 是待定参数
    lam = p[0]
    q = 10
    dy0 = y[1]
    dy1 = -(lam - 2 * q * np.cos(2 * x)) * y[0]
    return np.vstack((dy0, dy1))


# 计算 边界条件
def boundCond(ya, yb, p):
    lam = p[0]
    return np.array([ya[0] - 1, ya[0], yb[0]])


xa, xb = 0, np.pi  # 边界点 (xa,xb)
xini = np.linspace(xa, xb, 11)  # 确定 x 的初值
xSol = np.linspace(xa, xb, 100)  # 输出的网格节点

for k in range(5):
    A = 0.75 * k
    y0ini = np.cos(8 * xini)  # 设置 y0 的初值
    y1ini = -A * np.sin(8 * xini)  # 设置 y1 的初值
    yini = np.vstack((y0ini, y1ini))  # 确定 y=[y0,y1] 的初值
    res = solve_bvp(dydx, boundCond, xini, yini, p=[10])  # 求解 BVP
    y0 = res.sol(xSol)[0]  # 网格节点处的 y 值
    y1 = res.sol(xSol)[1]  # 网格节点处的 y 值
    plt.plot(xSol, y0, '--')
    plt.plot(xSol, y1, '-', label='A = {:.2f}'.format(A))

plt.xlabel("hbu")
plt.ylabel("y")
plt.title("Characteristic function of Mathieu equation")
plt.axis([0, np.pi, -5, 5])
plt.legend(loc='best')
plt.text(2, -4, "hao-hbu", color='whitesmoke')
plt.show()

在这里插入图片描述

三、非线性规划

1.非线性规划问题的描述

在这里插入图片描述
在这里插入图片描述

2.scipy.optimize.brent() 求解单变量无约束优化问题

非线性规划最简单的形式是一维搜索，一维搜索的常用方法是函数逼近法和区间收缩法。

brent() 函数是 SciPy.optimize 模块中求解单变量无约束优化问题最小值的首选方法。这是牛顿法和二分法的混合方法，既能保证稳定性又能快速收敛。
在这里插入图片描述
optimize.brent() 的主要参数：
func: callable f(x,args) 　　目标函数 f ( x ) f(x)f(x)，以函数形式表示，可以通过 *args 传递参数
args: tuple　　可选项，以 f(x,*args) 的形式将可变参数 p 传递给目标函数 f ( x , p ) f(x,p)f(x,p) 。
brack: tuple　　可选项，搜索算法的开始区间（不是指 x 的上下限）

optimize.brent() 的主要返回值：
xmin: 　　返回函数达到最小值时的 x（注意是局部最优，不一定是全局最优）。
fval: 　　返回函数的最优值（默认不返回，仅当 full_output 为 1 时返回）。

optimize.brent() 的使用例程：

import numpy as np
from scipy.optimize import brent
import matplotlib.pyplot as plt


# 1. Demo1：单变量无约束优化问题(Scipy.optimize.brent)
def objf(x):  # 目标函数
    fx = x ** 2 - 8 * np.sin(2 * x + np.pi)
    return fx


xIni = -5.0
xOpt = brent(objf, brack=(xIni, 2))
print("xIni={:.4f}\tfxIni={:.4f}".format(xIni, objf(xIni)))
print("xOpt={:.4f}\tfxOpt={:.4f}".format(xOpt, objf(xOpt)))

x = np.linspace(-10, 10, 500)
plt.plot(x, objf(x))
plt.xlabel("x")
plt.ylabel("objf(x)")
plt.text(xIni, objf(xIni), "Initial")
plt.text(xOpt, objf(xOpt), "Best")
plt.scatter(xIni, objf(xIni))
plt.scatter(xOpt, objf(xOpt))
plt.show()

运行结果：
在这里插入图片描述
在这里插入图片描述

3.scipy.optimize.fmin() 求解多变量无约束优化问题

fmin() 函数是 SciPy.optimize 模块中求解多变量无约束优化问题（最小值）的首选方法，采用下山单纯性方法。下山单纯性方法又称 Nelder-Mead 法，只使用目标函数值，不需要导数或二阶导数值，是最重要的多维无约束优化问题数值方法之一。
在这里插入图片描述
在这里插入图片描述

from scipy.optimize import brent, fmin, minimize
import numpy as np
​
# 2. Demo2：多变量无约束优化问题(Scipy.optimize.brent)
# Rosenbrock 测试函数
def objf2(x):  # Rosenbrock benchmark function
    fx = sum(100.0 * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0)
    return fx
​
xIni = np.array([-2, -2])
xOpt = fmin(objf2, xIni)
print("xIni={:.4f},{:.4f}\tfxIni={:.4f}".format(xIni[0],xIni[1],objf2(xIni)))
print("xOpt={:.4f},{:.4f}\tfxOpt={:.4f}".format(xOpt[0],xOpt[1],objf2(xOpt)))

在这里插入图片描述

4.scipy.optimize.minimize() 求解非线性规划问题

minimize() 函数是 SciPy.optimize 模块中求解多变量优化问题的通用方法，可以调用多种算法，支持约束优化和无约束优化。
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

import numpy as np
from scipy.optimize import minimize


# 定义目标函数
def objf3(x):  # Rosenbrock 测试函数
    fx = sum(100.0 * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0)
    return fx


# 定义边界约束（优化变量的上下限）
b0 = (0.0, None)  # 0.0 <= x[0] <= Inf
b1 = (0.0, 10.0)  # 0.0 <= x[1] <= 10.0
b2 = (-5.0, 100.)  # -5.0 <= x[2] <= 100.0
bnds = (b0, b1, b2)  # 边界约束

# 优化计算
xIni = np.array([1., 2., 3.])
resRosen = minimize(objf3, xIni, method='SLSQP', bounds=bnds)
xOpt = resRosen.x

print("xOpt = {:.4f}, {:.4f}, {:.4f}".format(xOpt[0], xOpt[1], xOpt[2]))
print("min f(x) = {:.4f}".format(objf3(xOpt)))

在这里插入图片描述

5.约束非线性规划问题实例

在这里插入图片描述
在这里插入图片描述

import numpy as np
from scipy.optimize import minimize


def objF4(x):  # 定义目标函数
    a, b, c, d = 1, 2, 3, 8
    fx = a * x[0] ** 2 + b * x[1] ** 2 + c * x[2] ** 2 + d
    return fx


# 定义约束条件函数
def constraint1(x):  # 不等式约束 f(x)>=0
    return x[0] ** 2 - x[1] + x[2] ** 2


def constraint2(x):  # 不等式约束 转换为标准形式
    return -(x[0] + x[1] ** 2 + x[2] ** 3 - 20)


def constraint3(x):  # 等式约束
    return -x[0] - x[1] ** 2 + 2


def constraint4(x):  # 等式约束
    return x[1] + 2 * x[2] ** 2 - 3


# 定义边界约束
b = (0.0, None)
bnds = (b, b, b)

# 定义约束条件
con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'ineq', 'fun': constraint2}
con3 = {'type': 'eq', 'fun': constraint3}
con4 = {'type': 'eq', 'fun': constraint4}
cons = ([con1, con2, con3, con4])  # 3个约束条件

# 求解优化问题
x0 = np.array([1., 2., 3.])  # 定义搜索的初值
res = minimize(objF4, x0, method='SLSQP', bounds=bnds, constraints=cons)

print("Optimization problem (res):\t{}".format(res.message))  # 优化是否成功
print("xOpt = {}".format(res.x))  # 自变量的优化值
print("min f(x) = {:.4f}".format(res.fun))  # 目标函数的优化值

在这里插入图片描述

相关文章：

  • 你把 《时间》 玩明白
  • 飞桨机器学习最小实现，出租车计费规则预测
  • 【andriod】设备APP开发之各种细节部署和操作
  • 矩阵论习题1.1
  • OAuth2-单点-多点-三方登录
  • SSO 基于token vue + element ui spring boot前端分离
  • pandas：jupyter notebook笔记（更新中）
  • 牛客网刷题——JAVA
  • Tpflow V6.0.8 正式版发布
  • 中缀转后缀表达式(思路分析) [数据结构][Java]
  • WhatsApp Business 与 Google My Business：它们的运作方式和不同之处
  • 【网络】网络基础知识点(MAC、IP、通讯模式、通信方式、集线器、交换机、路由器、冲突域/广播域、ARP、VLAN等)
  • Windows 应用商店无法打开解决办法
  • Jupyter Notebook与Pycharm代码在Docker容器中运行
  • 【组成原理-IO系统】IO控制方式
  • 【347天】每日项目总结系列085（2018.01.18）
  • CSS 三角实现
  • httpie使用详解
  • IDEA 插件开发入门教程
  • LeetCode18.四数之和 JavaScript
  • Linux学习笔记6-使用fdisk进行磁盘管理
  • NSTimer学习笔记
  • 第三十一到第三十三天：我是精明的小卖家（一）
  • 机器人定位导航技术 激光SLAM与视觉SLAM谁更胜一筹？
  • 如何优雅的使用vue+Dcloud（Hbuild）开发混合app
  • 使用 Docker 部署 Spring Boot项目
  • 算法-插入排序
  • 验证码识别技术——15分钟带你突破各种复杂不定长验证码
  • 在electron中实现跨域请求，无需更改服务器端设置
  • ​第20课 在Android Native开发中加入新的C++类
  • # 20155222 2016-2017-2 《Java程序设计》第5周学习总结
  • $.each()与$(selector).each()
  • （145）光线追踪距离场柔和阴影
  • (C++20) consteval立即函数
  • （C语言版）链表（三）——实现双向链表创建、删除、插入、释放内存等简单操作...
  • （Redis使用系列） Springboot 整合Redisson 实现分布式锁 七
  • （SpringBoot）第七章：SpringBoot日志文件
  • （十三）Flask之特殊装饰器详解
  • (原創) 博客園正式支援VHDL語法著色功能 (SOC) (VHDL)
  • （转）JVM内存分配 -Xms128m -Xmx512m -XX:PermSize=128m -XX:MaxPermSize=512m
  • (转)拼包函数及网络封包的异常处理(含代码)
  • ＊Django中的Ajax 纯js的书写样式1
  • .NET 2.0中新增的一些TryGet,TryParse等方法
  • .NET 编写一个可以异步等待循环中任何一个部分的 Awaiter
  • .sdf和.msp文件读取
  • @DataRedisTest测试redis从未如此丝滑
  • @Documented注解的作用
  • @Pointcut 使用
  • @RequestMapping用法详解
  • [ 2222 ]http://e.eqxiu.com/s/wJMf15Ku
  • [ C++ ] STL---string类的模拟实现
  • [ vulhub漏洞复现篇 ] ECShop 2.x / 3.x SQL注入/远程执行代码漏洞 xianzhi-2017-02-82239600
  • [ 云计算 | Azure 实践 ] 在 Azure 门户中创建 VM 虚拟机并进行验证
  • [3300万人的聊天室] 作为产品的上游公司该如何？
  • [Android]如何调试Native memory crash issue