May 28, 2026

miniyuan

函数逼近

---

题 6.1

取倒数得：

$$\frac{1}{y}=a+\frac{b}{x}$$

令 $Y=\dfrac{1}{y},\;X=\dfrac{1}{x}$，则化为线性模型 $Y=a+bX$。

列出正规方程组：

$$\begin{pmatrix} 5 & \sum_{i=1}^5 X_i \\ \sum_{i=1}^5 X_i & \sum_{i=1}^5 X_i^2 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} \sum_{i=1}^5 Y_i \\ \sum_{i=1}^5 X_i Y_i \end{pmatrix}$$

也即：

$$\begin{pmatrix} 5 & 2.2833 \\ 2.2833 & 1.4636 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} 12.2841 \\ 6.0326 \end{pmatrix}$$

解得：

$$a \approx 1.9980,\quad b \approx 1.0048$$

故拟合函数为：

$$y=\frac{x}{1.9980x+1.0048}$$

---

题 6.2

令：

$$A=\begin{pmatrix}1&-1\\-1&1\\2&-2\end{pmatrix},\qquad \boldsymbol{b}=\begin{pmatrix}5\\-4\\10\end{pmatrix}$$

则正规方程为：

$$A^{\mathsf T}A\boldsymbol{x} = A^{\mathsf T}\boldsymbol{b}$$

其中：

$$A^{\mathsf T}A=\begin{pmatrix}6&-6\\-6&6\end{pmatrix},\qquad A^{\mathsf T}\boldsymbol{b}=\begin{pmatrix}29\\-29\end{pmatrix}$$

注意矩阵 $A^{\mathsf T}A$ 奇异。从而最小二乘解不唯一，只能解得如下形式通解：

$$x_1=x_2+\frac{29}{6},\qquad x_2\in\mathbb{R}$$

考虑取其中的最小范数解，也即：

$$\min (x_1^2+x_2^2) = \min [2x_2^2+\frac{29}{3}x_2+(\frac{29}{6})^2]$$

解得最小范数解：

$$x_1=\frac{29}{12}\approx 2.4167,\qquad x_2=-\frac{29}{12}\approx -2.4167$$

---

题 6.3

取基函数 $\varphi_1(x)=\dfrac{1}{x},\;\varphi_2(x)=x$，定义离散内积 $(f,g)=\sum_{i=1}^{5}f(x_i)g(x_i)$。

则有：

$$\begin{aligned} (\varphi_1,\varphi_1)&=\sum\frac{1}{x_i^2}=1.3559 \\ (\varphi_1,\varphi_2)&=\sum 1=5 \\ (\varphi_2,\varphi_2)&=\sum x_i^2=121 \\ (y,\varphi_1)&=\sum\frac{y_i}{x_i}=16.0654 \\ (y,\varphi_2)&=\sum x_iy_i=304.5100 \end{aligned}$$

列出正规方程组：

$$\begin{pmatrix} 1.3559 & 5 \\ 5 & 121 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} 16.0654 \\ 304.5100 \end{pmatrix}$$

解得：

$$a\approx 3.0300,\qquad b\approx 2.3914$$

从而最佳平方逼近为：

$$y=\frac{3.0300}{x}+2.3914x$$

---

题 6.4

前三阶 Legendre 正交多项式为：

$$P_0(x)=1,\qquad P_1(x)=x,\qquad P_2(x)=\frac{3x^2-1}{2}$$

且内积 $(P_k,P_k)=\dfrac{2}{2k+1}$。

设最佳平方逼近多项式为 $S(x)=c_0P_0(x)+c_1P_1(x)+c_2P_2(x)$，$f(x)=\cos(\pi x)$，则展开系数：

$$c_k=\frac{(f,P_k)}{(P_k,P_k)} =\frac{2k+1}{2}(f,P_k) =\frac{2k+1}{2}\int_{-1}^{1}\cos(\pi x)P_k(x)\,\mathrm{d}x$$

从而：

$$c_0 =\frac{1}{2}\int_{-1}^{1}\cos(\pi x)\,\mathrm{d}x =\frac{\sin(\pi x)}{2\pi}\Big|_{-1}^{1} =0$$ $$c_1 =\frac{3}{2}\int_{-1}^{1}\cos(\pi x)x\,\mathrm{d}x =0$$ $$c_2 =\frac{5}{2}\int_{-1}^{1}\cos(\pi x)\frac{3x^2-1}{2}\,\mathrm{d}x =\frac{15x\cos(\pi x)}{2\pi^2}\Big|_{-1}^{1} =-\frac{15}{\pi^2}$$

从而最佳平方逼近多项式为：

$$\begin{aligned} S(x)&=c_0P_0(x)+c_1P_1(x)+c_2P_2(x) \\ &=-\frac{15}{\pi^2}\cdot\frac{3x^2-1}{2} \\ &=-\frac{45}{2\pi^2}x^2+\frac{15}{2\pi^2} \\ &\approx -2.2797\,x^2+0.7599 \end{aligned}$$

题 6.5

$$\cos^n \theta = (\frac{\mathrm{e}^{i\theta} + \mathrm{e}^{-i\theta}}{2})^n = \frac{1}{2^n} \sum_{k=0}^n \mathrm{C}_n^k \mathrm{e}^{i(n-2k)\theta}$$

且有：

$$\cos^n \theta = \cos^n(-\theta) = \frac{1}{2^n} \sum_{k=0}^n \mathrm{C}_n^k \mathrm{e}^{-i(n-2k)\theta}$$

相加即得：

$$\cos^n \theta = \frac{1}{2^n} \sum_{k=0}^n \mathrm{C}_n^k \frac{\mathrm{e}^{i(n-2k)\theta} + \mathrm{e}^{-i(n-2k)\theta}}{2} = \frac{1}{2^n} \sum_{k=0}^n \mathrm{C}_n^k \cos((n-2k)\theta)$$

令 $x=\cos\theta$，由于 $T_n(x)=T_n(\cos\theta)=\cos(n\theta)$ 代入即得：

$$x^n = \frac{1}{2^n} \sum_{k=0}^n \mathrm{C}_n^k \cos(|n-2k|\theta) = \frac{1}{2^n} \sum_{k=0}^n \mathrm{C}_n^k T_{|n-2k|}(x)$$

从而当 $n$ 为奇数时：

$$2^{n-1}c_{n-2j} = \mathrm{C}_n^j,\quad j=0,1,\ldots,\frac{n-1}{2}$$

当 $n$ 为偶数时：

$$\begin{aligned} &2^{n-1}c_{n-2j} = \mathrm{C}_n^j,\quad j=0,1,\ldots,\frac{n}{2}-1 \\[6pt] &2^{n-1}c_0 = \frac{1}{2}\mathrm{C}_n^{n/2} \end{aligned}$$

运行结果：

input: 6
output: 1 6 15 10

input: 7
output: 1 7 21 35

input: 11
output: 1 11 55 165 330 462

input: 12
output: 1 12 66 220 495 792 462

input: 15
output: 1 15 105 455 1365 3003 5005 6435

题 6.6

使用 SINDy 进行非线性动力学的稀疏辨识。

1. 数据标准化

先对数据做标准化，避免不同族群间数量级差异问题。对每个物种 $i$ 的观测序列：

$$\begin{aligned} \mu_i&=\frac{1}{m+1}\sum_{j=0}^{m}p_{i,j} \\ \sigma_i&=\sqrt{\frac{1}{m}\sum_{j=0}^{m}(p_{i,j}-\mu_i)^2} \\ \tilde p_{i,j}&=\frac{p_{i,j}-\mu_i}{\sigma_i} \end{aligned}$$

后续所有运算均在标准化后的数据 $\{\tilde p_{i,j}\}$ 上进行。

若需还原到原量纲，可利用 $\tilde p_i=(p_i-\mu_i)/\sigma_i$ 做变量替换，将方程改写为关于原始种群数量 $p_i$ 的表达式。

2. 导数数据构造

对每个物种 $i$ 和时刻 $j=0,1,\dots,m-1$：

$$\dot{\tilde p}_{i,j}=\frac{\tilde p_{i,j+1}-\tilde p_{i,j}}{\Delta t}$$

得到 $m$ 个样本点。

3. 构建候选函数库

假设 $f_i$ 可以表示为以下多元多项式基的线性组合：

$$\Theta(\boldsymbol p)=\begin{bmatrix} 1,\; p_1,\; p_2,\; \dots,\; p_n,\; p_1^2,\; p_1p_2,\; \dots,\; p_n^2,\; \dots \end{bmatrix}^T$$

不妨取多项式最高阶数为 $d=2$，则库中共有 $K=\dbinom{n+d}{d}=\dfrac{(n+2)(n+1)}{2}$ 个基函数，记为：

$$\theta_1(\boldsymbol p),\;\theta_2(\boldsymbol p),\;\dots,\;\theta_K(\boldsymbol p)$$

其中线性项可反映种群自然增长/衰减；二次项 $p_kp_l$ 可反映种间相互作用（捕食、竞争、互惠）；常数项可反映外部输入。

4. 建立回归模型

对每个物种 $i$，假设：

$$f_i(\boldsymbol p)=\sum_{k=1}^{K}\xi_{i,k}\,\theta_k(\boldsymbol p)$$

其中 $\xi_{i,k}$ 为待求系数。代入所有 $m$ 个采样时刻，得到线性方程组：

$$\underbrace{\begin{bmatrix} \dot{\tilde p}_{i,0}\\ \dot{\tilde p}_{i,1}\\ \vdots\\ \dot{\tilde p}_{i,m-1} \end{bmatrix}}_{\displaystyle \dot{\tilde{\mathbf P}}_i\in\mathbb R^{m}} = \underbrace{\begin{bmatrix} \theta_1(\tilde{\boldsymbol p}_0) & \theta_2(\tilde{\boldsymbol p}_0) & \cdots & \theta_K(\tilde{\boldsymbol p}_0)\\ \theta_1(\tilde{\boldsymbol p}_1) & \theta_2(\tilde{\boldsymbol p}_1) & \cdots & \theta_K(\tilde{\boldsymbol p}_1)\\ \vdots & \vdots & \ddots & \vdots\\ \theta_1(\tilde{\boldsymbol p}_{m-1}) & \theta_2(\tilde{\boldsymbol p}_{m-1}) & \cdots & \theta_K(\tilde{\boldsymbol p}_{m-1}) \end{bmatrix}}_{\displaystyle \boldsymbol\Theta(\tilde{\mathbf P})\in\mathbb R^{m\times K}} \underbrace{\begin{bmatrix} \xi_{i,1}\\ \xi_{i,2}\\ \vdots\\ \xi_{i,K} \end{bmatrix}}_{\displaystyle \boldsymbol\xi_i\in\mathbb R^{K}}$$

简记为：

$$\dot{\tilde{\mathbf P}}_i = \boldsymbol\Theta(\tilde{\mathbf P})\,\boldsymbol\xi_i$$

5. 最小二乘求解

由于 $m\gg K$，这是超定线性系统，可用最小二乘法求解：

$$\boldsymbol\xi_i=\arg\min_{\boldsymbol\xi}\bigl\|\dot{\tilde{\mathbf P}}_i-\boldsymbol\Theta(\tilde{\mathbf P})\boldsymbol\xi\bigr\|_2^2$$

正规方程组解为：

$$\boldsymbol\xi_i =\bigl(\boldsymbol\Theta(\tilde{\mathbf P})^\mathsf T\boldsymbol\Theta(\tilde{\mathbf P})\bigr)^{-1} \boldsymbol\Theta(\tilde{\mathbf P})^\mathsf T\dot{\tilde{\mathbf P}}_i$$

6. 稀疏正则化

在实际的生态系统中，每个物种通常只与少数其他物种存在直接相互作用。 因此 $\boldsymbol\xi_i$ 应该是稀疏向量。

最小二乘解通常非稀疏且过拟合，从而可以改用 LASSO：

$$\min_{\boldsymbol\xi_i}\bigl\|\dot{\tilde{\mathbf P}}_i-\boldsymbol\Theta\boldsymbol\xi_i\bigr\|_2^2+\lambda\|\boldsymbol\xi_i\|_1$$

其中 $\lambda>0$ 为正则化超参数，控制稀疏度。

$\lambda$ 可通过 K-折交叉验证选取：在训练集上求解 $\boldsymbol\xi_i$，在验证集上评估预测误差，选择使验证误差最小的 $\lambda$。

7. 组装微分方程

对第 $i$ 个物种，设经 LASSO 求解后得到的稀疏系数为 $\hat\xi_{i,1},\dots,\hat\xi_{i,K}$。 则该物种的演化方程为：

$$\frac{\mathrm{d}\tilde p_i}{\mathrm{d}t} =f_i(\tilde p_1,\dots,\tilde p_n)=\sum_{k=1}^K \hat\xi_{i,k}\,\theta_k(\tilde{\boldsymbol p}) ,\qquad i=1,2,\dots,n$$

将 $n$ 个物种的方程联立，即得到完整的生态系统常微分方程组。