DW数学建模学习笔记—

第1章解析方法与几何模型

1.1 向量表示法与几何建模基本案例

1.1.1 几何建模的思想

数学世界的探索始于计数和丈量。计数带来了数字和代数学，而丈量则推动了几何学的发展。高中毕业后，我们对数学模型的理解还较浅，但几何模型对我们而言是最直观的。通过几何定理，我们可以直观地感受几何关系，并利用这些关系推算长度等。

几何模型可以分为位置关系和数量关系。位置关系包括平行、垂直等，数量关系涉及边长、角度、面积等求解。分析几何问题的方法主要有三种：

传统几何的演绎-证明体系：基于已证明的公理和定理，通过逻辑推理得出结果。
基于向量的计算化几何：将几何问题转化为向量的计算问题，方便求解。
基于极坐标与方程的解析几何：利用极坐标和参数方程，将几何问题转化为代数问题求解。

常见的几何定理包括三角形内角和定理、勾股定理、正弦定理、余弦定理、圆幂定理等。这些定理在解决几何问题时非常有用。

1.1.2 向量表示与坐标变换

向量是几何中表示方向和距离的工具，可以在高维空间中进行加减运算、数量乘运算等。在Python中，可以使用NumPy库创建和操作向量。例如：

import numpy as np
x = np.array([1, 2, 3, 5, 8])

向量的引入不仅仅是为了表示几何图形中的方向和距离，更重要的是利用代数方法解决几何问题。在物理问题、计算机图形学等领域，向量运算简化了计算过程。

二维空间的旋转变换：

import numpy as nptheta = np.radians(30)
rotation_matrix = np.array([[np.cos(theta), -np.sin(theta)],[np.sin(theta), np.cos(theta)]
])
point = np.array([a, b])
rotated_point = rotation_matrix.dot(point)
print("旋转后的坐标为:", rotated_point)

三维空间的旋转变换：

import numpy as npalpha = np.radians(30)
beta = np.radians(45)
gamma = np.radians(60)R_z = np.array([[np.cos(alpha), -np.sin(alpha), 0],[np.sin(alpha), np.cos(alpha), 0],[0, 0, 1]])R_y = np.array([[np.cos(beta), 0, np.sin(beta)],[0, 1, 0],[-np.sin(beta), 0, np.cos(beta)]])R_x = np.array([[1, 0, 0],[0, np.cos(gamma), -np.sin(gamma)],[0, np.sin(gamma), np.cos(gamma)]])R = R_z @ R_y @ R_xP = np.array([1, 2, 3])
P_rotated = R @ Pprint("旋转后P点的坐标为:", P_rotated)

1.2 Numpy 与线性代数

1.2.1 Numpy向量与矩阵的操作

在科学计算中，Numpy的数组对象是处理问题的得力助手。它能进行向量化运算，提升代码简洁性和运行速度。

创建向量和矩阵：

import numpy as np
vector = np.array([1, 2, 3])
matrix = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

向量和矩阵的基本属性：

print(vector.shape)  # (3,)
print(matrix.shape)  # (3, 3)

索引和切片：

print(vector[0])  # 输出第一个元素, 1
print(matrix[1, 1])  # 输出第二行第二列的元素, 5

向量和矩阵的运算：

vector1 = np.array([1, 2, 3])
vector2 = np.array([4, 5, 6])
print(np.add(vector1, vector2))  # [5, 7, 9]matrix1 = np.array([[1, 2], [3, 4]])
matrix2 = np.array([[5, 6], [7, 8]])
print(np.dot(matrix1, matrix2))  # [[19, 22], [43, 50]]

1.2.2 利用Numpy进行线性代数基本运算

Numpy的广播机制和线性代数函数使得数值计算更加便捷。

数量乘法：

scalar = 5
scaled_vector = scalar * vector
print("Scaled vector:", scaled_vector)

矩阵的转置：

transposed_matrix = matrix.T
print("Transposed matrix:\n", transposed_matrix)

行列式：

matrix_determinant = np.linalg.det(matrix)
print("Matrix determinant:", matrix_determinant)

求解线性方程组：

A = np.array([[3, 1], [1, 2]])
b = np.array([9, 8])
solution = np.linalg.solve(A, b)
print("Solution of the linear system:", solution)

1.2.3 numpy.linalg 的使用

Numpy的 linalg 子模块提供了丰富的线性代数函数。

计算逆矩阵：

pseudo_inverse_matrix = np.linalg.pinv(matrix)
print("Pseudo-inverse of the matrix:\n", pseudo_inverse_matrix)

特征值和特征向量：

eigenvalues, eigenvectors = np.linalg.eig(matrix)
print(eigenvalues)
print(eigenvectors)

奇异值分解：

U, S, V = np.linalg.svd(matrix)
print(U)
print(S)
print(V)

范数计算：

norm = np.linalg.norm(vector)
print(norm)

1.3 平面几何模型的构建

1.3.1 问题背景介绍

以2023年高教社杯全国大学生数学建模竞赛B题为例，单波束和多波束测深系统用于测量海底深度。多波束系统在垂直于航迹的平面内发射多个波束，覆盖更大面积。我们需要建立多波束测深的覆盖宽度及相邻条带之间重叠率的数学模型。

1.3.2 问题一的分析

问题一是几何模型问题，通过三角函数和几何关系推导出覆盖宽度和重叠率。

1.3.3 问题一的模型建立

构造几何模型，如图所示，利用三角形中的正弦定理和几何关系，推导出覆盖宽度和重叠部分的计算方法。

几何模型推导：

import numpy as np
from scipy.optimize import fsolvetheta = 2 * np.pi / 3
alpha = 1.5 / 180 * np.pi
htheta = theta / 2
h = 70
d = 200
k = np.tan(np.pi / 2 - htheta)
k0 = np.tan(alpha)Aleft = []
Aright = []
Acenter = []
W = []for n in range(-4, 5):leftsolve = lambda t: k * (t - n * d) - k0 * t + hrightsolve = lambda t: -k * (t - n * d) - k0 * t + htleft = fsolve(leftsolve, 0)tright = fsolve(rightsolve, 0)Aleft.append([tleft[0], k0 * tleft[0] - h])Aright.append([tright[0], k0 * tright[0] - h])Acenter.append([200 * n, k0 * 200 * n - h])
Aleft = np.array(Aleft)
Aright = np.array(Aright)
Acenter = np.array(Acenter)
D = Acenter[:, 1]
W = np.sqrt((Aleft[:, 0] - Aright[:, 0]) ** 2 + (Aleft[:, 1] - Aright[:, 1]) ** 2)cover = np.zeros(8)
for i in range(8):cover[i] = np.sqrt((Aright[i, 0] - Aleft[i + 1, 0]) ** 2 + (```pythonAright[i, 1] - Aleft[i + 1, 1]) ** 2)
eta = cover / W[1:]print("海水深度 D:", D)
print("覆盖宽度 W:", W)
print("重合部分比例 eta:", eta)

1.3.4 问题一的模型求解与讨论

通过建模与编程求解，可以得到如下计算结果：

$d/\text{m}$	$- 800$	$- 600$	$- 400$	$- 200$	$0$	$200$	$400$	$600$	$800$
$D/\mathrm{m}$	$90.94$	$85.711$	$80.474$	$75.237$	$70$	$64.76$	$59.526$	$54.289$	$49.051$
$W/\mathrm{m}$	$315.81$	$297.62$	$279.441$	$261.25$	$243$	$224.8$	$206.698$	$188.513$	$170.328$
$\eta/\%$	—	$0.3569$	$0.3151$	$0.267$	$0.21$	$0.148$	$0.0740$	$0.0153$	$- 0.12365$

从表中可以看出，测线之间的重叠率在不同距离下有所变化。通过Python编程可以更高效地进行这种几何计算。

1.4 立体几何模型的构建

1.4.1 问题背景介绍

继续以2023年高教社杯全国大学生数学建模竞赛B题为例，考虑一个矩形待测海域，测线方向与海底坡面的法向在水平面上投影的夹角为 $\beta$ 。我们需要建立多波束测深覆盖宽度的数学模型。

1.4.2 问题二的分析

问题二是问题一的延续，将场景从二维向三维扩展，通过构造四棱锥几何模型，利用三面角定理求解。

1.4.3 问题二的模型建立

构造三维几何模型，通过几何关系和三面角公式，推导覆盖宽度的计算方法。

import numpy as np
from scipy.optimize import fsolvetheta = 2 * np.pi / 3
alpha = -1.5 / 180 * np.pi
htheta = theta / 2
h = 120
unit = 1852
k0 = np.tan(alpha)W = np.zeros((8, 8))for i in range(1, 9):for j in range(1, 9):beta = (i - 1) * np.pi / 4d = (j - 1) * 0.3 * unitv = np.array([np.cos(beta), np.sin(beta), 0])origin = v * dv1 = np.array([-np.sin(beta) * np.sin(htheta), np.cos(beta) * np.sin(htheta), -np.cos(htheta)])v2 = np.array([np.sin(beta) * np.sin(htheta), -np.cos(beta) * np.sin(htheta), -np.cos(htheta)])leftsolve = lambda t: (v1[0] * t + origin[0]) * k0 - h - (v1[2] * t + origin[2])rightsolve = lambda t: (v2[0] * t + origin[0]) * k0 - h - (v2[2] * t + origin[2])tleft = fsolve(leftsolve, 0)tright = fsolve(rightsolve, 0)pleft = v1 * tleft + originpright = v2 * tright + originW[i-1, j-1] = np.linalg.norm(pleft - pright)print("覆盖宽度矩阵 W:\n", W)

1.4.4 问题二的模型求解与讨论

通过建模与编程求解，可以得到如下计算结果：

测量船距海域中心点处的距离 (m)	覆盖宽度 (m)	测线方向夹角 (度)
0	415.69	0
0.3 * 1852	466.09	45
0.6 * 1852	516.49	90
0.9 * 1852	566.89	135
1.2 * 1852	617.29	180
1.5 * 1852	667.69	225
1.8 * 1852	718.09	270
2.1 * 1852	768.49	315

从表中可以看出，测量船距海域中心点的距离与覆盖宽度之间存在一定关系，随着距离的增加，覆盖宽度逐渐增加。

1.5 使用Python解方程与方程组

1.5.1 利用Numpy求线性方程（组）的数值解

Numpy库提供了强大的线性代数运算功能，可以方便地求解线性方程组。

import numpy as npa = np.array([[10, -1, -2], [-1, 10, -2], [-1, -1, 5]])
b = np.array([[72], [83], [42]])
c = np.linalg.solve(a, b)
print(c)# 或使用矩阵的逆来求解
x = np.linalg.inv(a).dot(b)
print(x)

1.5.2 利用Sympy求方程（组）的解析解

Sympy库提供了符号计算功能，可以求解方程和方程组的解析解。

from sympy import symbols, solve, nonlinsolvex, y = symbols('x y')
print(solve(x * 2 - 2, x))  # 解方程2x - 2 = 0
print(solve([x + y - 35, x * 2 + y * 4 - 94], x, y))  # 解方程组x + y = 35, 2x + 4y = 94
print(solve(x**2 + x - 20, x))  # 解方程x^2 + x - 20 = 0a, b, c, d = symbols('a b c d', real=True)
print(nonlinsolve([a**2 + a + b, a - b], [a, b]))  # 解非线性方程组a^2 + a + b = 0, a - b = 0

1.5.3 利用Scipy求方程（组）的数值解

对于复杂的非线性方程组，可以使用Scipy库中的fsolve函数求解。

from scipy.optimize import fsolve
from math import sin, cos, pidef equations(vars):x, y, theta = varsL1, L2, L3 = 3, 3, 3p1, p2, p3 = 5, 5, 3x1, x2, y2 = 5, 0, 6eq1 = (x + L3*cos(theta) - x1)**2 + (y + L3*sin(theta))**2 - p2**2eq2 = x**2 + y**2 - p1**2eq3 = (x + L2*cos(pi/3 + theta))**2 + (y + L2*sin(pi/3 + theta) - y2)**2 - p3**2return [eq1, eq2, eq3]initial_guess = [-1.37, 4.80, 0.12]
result = fsolve(equations, initial_guess)
print(result)

第2章微分方程与动力系统

2.1 微分方程的理论基础

微分方程在数学和工程中有着广泛的应用。微分方程描述了一个函数及其导数之间的关系，通过这种关系我们可以研究和求解许多实际问题。

2.1.1 函数、导数与微分

微分方程的核心在于微分和导数的概念。导数描述了函数在某一点的变化率，微分则描述了当自变量施加微小增量时函数值的变化。当增量非常小时，函数的变化量接近于该点处切线的变化量。

公式 (2.1.1) 展示了微分的基本形式：
[ \frac{d y}{d x} = f’(x) ]

2.1.2 一阶线性微分方程的解

一阶线性微分方程的一般形式为：
[ \frac{d y}{d x} + y P(x) = Q(x) ]

通过分离变量法和常数变易法可以求得其解。对于齐次方程：
[ \frac{d y}{d x} + y P(x) = 0 ]

其解的通式为：
[ y = C \exp \left( -\int P(x) , dx \right) ]

对于非齐次方程的通解则为：
[ y = \exp \left( -\int P(x) , dx \right) \left[ \int Q(x) \exp \left( \int P(x) , dx \right) , dx + C \right] ]

2.1.3 二阶常系数线性微分方程的解

二阶常系数线性微分方程的一般形式为：
[ f’‘(x) + p f’(x) + q f(x) = C(x) ]

通过特征根法求解齐次方程：
[ r^2 + pr + q = 0 ]

特征根的类型决定了解的形式。对于非齐次方程，可以通过特解与齐次方程通解的组合得到总解。

2.1.4 利用Python求函数的微分与积分

使用Python的Numpy和SciPy库可以方便地计算函数的微分和积分。例如，SciPy库中的quad函数可以求定积分：

import numpy as np
from scipy.integrate import quaddef f(x):return np.cos(2 * np.pi * x) * np.exp(-x) + 1.2integral, error = quad(f, 0, 0.7)
print(f'定积分的结果是：{integral}')

通过梯形法则近似计算函数的定积分：

h = x[1] - x[0]
xn = 0.7
s = 0
for i in range(1000):xn1 = xn + hyn = np.cos(2 * np.pi * xn) * np.exp(-xn) + 1.2yn1 = np.cos(2 * np.pi * xn1) * np.exp(-xn1) + 1.2s0 = (yn + yn1) * h / 2s += s0xn = xn1
s

2.2 使用SciPy和Sympy解微分方程

多数微分方程没有解析解，因此使用数值方法求解是常见的做法。

2.2.1 使用sympy求解微分方程解析解

使用sympy库中的dsolve函数可以求解微分方程的解析解。例如：

from sympy import symbols, Function, Eq, dsolvey = symbols('y', cls=Function)
x = symbols('x')
eq = Eq(y(x).diff(x, 2) + 2*y(x).diff(x, 1) + y(x), x**2)
dsolve(eq, y(x))

2.2.2 使用scipy求解微分方程数值解

使用scipy库中的odeint和solve_ivp函数可以求解微分方程的数值解。例如：

import numpy as np
from scipy.integrate import odeintdy = lambda y, x: 1/(1 + x**2) - 2*y**2
x = np.arange(0, 10.5, 0.1)
sol = odeint(dy, 0, x)
print("x={}\n对应的数值解y={}".format(x, sol.T))

2.3 偏微分方程的数值求解

偏微分方程用于描述多元函数及其偏导数之间的关系。通过离散化方法，如有限差分法，可以求解偏微分方程。

2.3.1 偏微分方程数值解的理论基础

偏微分方程可以通过离散化转换为代数方程组，然后通过数值方法求解。常用的方法包括有限差分方法、有限元方法等。

2.3.2 偏微分方程数值解的应用案例

应用有限差分法求解热传导方程和波动方程等实际问题。例如：

import numpy as np
import matplotlib.pyplot as plth = 0.1  # 空间步长
N = 30  # 空间步数
dt = 0.0001  # 时间步长
M = 10000  # 时间的步数
A = dt / (h ** 2)  # λτ/h^2
U = np.zeros([N+1, M+1])  # 建立二维空数组
Space = np.arange(0, (N+1)*h, h)  # 建立空间等差数列，从0到3，公差是h# 初始条件
for i in np.arange(0, N):U[i, 0] = 4 * i * h * (3 - i * h)# 递推关系
for k in np.arange(0, M):for i in np.arange(1, N):U[i, k+1] = A * U[i+1, k] + (1 - 2*A) * U[i, k] + A * U[i-1, k]plt.plot(Space, U[:, 0], 'g-', label='t=0')
plt.plot(Space, U[:, 3000], 'b-', label='t=3/10')
plt.plot(Space, U[:, 6000], 'k-', label='t=6/10')
plt.plot(Space, U[:, 9000], 'r-', label='t=9/10')
plt.plot(Space, U[:, 10000], 'y-', label='t=1')
plt.ylabel('u(x,t)')
plt.xlabel('x')
plt.legend(loc='upper right')
plt.show()

2.4 微分方程的应用案例

2.4.1 人口增长模型

人口增长可以用微分方程建模。Malthus模型和Logistic模型分别描述了简单的指数增长和受限增长。

Malthus模型：
[ \frac{d x}{d t} = r x, \quad x(t_0) = x_0 ]

Logistic模型：
[ \frac{d x}{d t} = r x \left(1 - \frac{x}{x_m}\right), \quad x(t_0) = x_0 ]

Python实现示例：

import pandas as pd
import numpy as np
from scipy import optimize
import matplotlib.pyplot as pltdata = pd.read_excel("总人口.xls")
y = data['年末总人口(万人)'].values.astype(float)
x_0 = y[0]def f2(x, r, K):return K / (1 + (K/x_0 - 1) * np.exp(-r * x))def f(x, r):return x_0 * np.exp(r * x)x = np.arange(0, 30, 1)
fita, _ = optimize.curve_fit(f, x, y)
plt.plot(np.arange(1990, 2020, 1), 10000 * y, '.')
plt.plot(np.arange(1990, 2020, 1), 10000 * f(x, fita[0]), 'r*-')
fita, _ = optimize.curve_fit(f2, x, y)
plt.plot(np.arange(1990, 2020, 1), 10000 * f2(x, fita[0], 11019), 'b^-')
plt.legend(["Origin", "Malthus", "Logistic"])
plt.xlabel("Year")
plt.ylabel("Population")
plt.title('Population Prediction')
plt.show()