miniyuan

数值微分与数值积分

题 7.1

精确值为：

f′(x)=−2(1+x)3f'(x)=-\frac{2}{(1+x)^3}

从而：

f′(0.5)=−1627≈−0.5925926f'(0.5)=-\frac{16}{27}\approx -0.5925926

取步长 h=0.1h=0.1 得函数表：

xix_i0.40.50.6
f(xi)f(x_i)2549\dfrac{25}{49}49\dfrac{4}{9}2564\dfrac{25}{64}

向前差商

f′(0.5)≈f(0.6)−f(0.5)0.1=−155288≈−0.5381944f'(0.5) \approx\frac{f(0.6)-f(0.5)}{0.1} =-\frac{155}{288} \approx-0.5381944

截断误差为：

e1(h)=−h2f′′(ξ)e_1(h) =-\frac{h}{2}f''(\xi)

其中 h=0.1h=0.1，ξ∈(0.5,0.6)\xi \in (0.5, 0.6)，f′′(x)=6(1+x)4f''(x)=\dfrac{6}{(1+x)^4}，故：

e1(h)=−310(1+ξ)4∈[−0.0593, −0.0458]e_1(h) =-\frac{3}{10(1+\xi)^4} \in[-0.0593,\,-0.0458]

向后差商

f′(0.5)≈f(0.5)−f(0.4)0.1=−290441≈−0.6575964f'(0.5) \approx\frac{f(0.5)-f(0.4)}{0.1} =-\frac{290}{441} \approx-0.6575964

截断误差为：

e2(h)=h2f′′(ξ)e_2(h) =\frac{h}{2}f''(\xi)

其中 h=0.1h=0.1，ξ∈(0.4,0.5)\xi \in (0.4,0.5)，f′′(x)=6(1+x)4f''(x)=\dfrac{6}{(1+x)^4}，故：

e2(h)=310(1+ξ)4∈[0.0593, 0.0781]e_2(h) =\frac{3}{10(1+\xi)^4} \in[0.0593,\,0.0781]

中心差商

f′(0.5)≈f(0.6)−f(0.4)0.2=−18753136≈−0.5978954f'(0.5) \approx\frac{f(0.6)-f(0.4)}{0.2} =-\frac{1875}{3136} \approx-0.5978954

截断误差为：

e3(h)=−h26f′′′(ξ)e_3(h) =-\frac{h^2}{6}f'''(\xi)

其中 h=0.1h=0.1，ξ∈(0.4,0.6)\xi \in (0.4,0.6)，f′′′(x)=−24(1+x)5f'''(x)=-\dfrac{24}{(1+x)^5}，故：

e3(h)=125(1+ξ)5∈[0.0038, 0.0074]e_3(h) =\frac{1}{25(1+\xi)^5} \in[0.0038,\,0.0074]

题 7.2

精确值为：

∫01e−x dx=(−e−x)01=1−e−1≈0.6321206\int_0^1 e^{-x}\,dx =(-e^{-x})_0^1 =1-e^{-1} \approx 0.6321206

梯形公式

∫01e−x dx≈1−02(e−0+e−1)=e+12e≈0.6839397\int_0^1 e^{-x}\,dx \approx \frac{1-0}{2}(e^{-0}+e^{-1}) =\frac{e+1}{2e} \approx 0.6839397

截断误差：

e1=−(1−0)312f′′(ξ)e_1 =-\frac{(1-0)^3}{12}f''(\xi)

其中 ξ∈[0,1]\xi \in [0, 1]，f′′(x)=e−xf''(x)=e^{-x}，故：

e1=−112e−ξ∈[−112, −112e]e_1 =-\frac{1}{12}e^{-\xi} \in[-\frac{1}{12},\,-\frac{1}{12e}]

也即 e1∈[−0.08333, −0.03066]e_1 \in [-0.08333,\,-0.03066]。

Simpson 公式

∫01e−x dx≈1−06(e−0+4e−0.5+e−1)=e+4e+16e≈0.6323337\int_0^1 e^{-x}\,dx \approx \frac{1-0}{6}(e^{-0}+4e^{-0.5}+e^{-1}) =\frac{e+4\sqrt{e}+1}{6e} \approx 0.6323337

截断误差：

e2=−(1−0)52880f(4)(ξ)e_2 =-\frac{(1-0)^5}{2880}f^{(4)}(\xi)

其中 ξ∈[0,1]\xi \in [0,1]，f(4)(x)=e−xf^{(4)}(x)=e^{-x}，故：

e2=−12880e−ξ∈[−12880, −12880e]e_2 =-\frac{1}{2880}e^{-\xi} \in [-\frac{1}{2880},\,-\frac{1}{2880e}]

也即 e2∈[−0.0003472, −0.0001277]e_2 \in [-0.0003472,\,-0.0001277]。

题 7.3

令求积公式对 f(x)=1, x, x2f(x)=1,\,x,\,x^2 成立，可得方程组：

{2h=A1+A2+A30=−h2A1+h2A32h33=h24A1+h24A3\begin{cases} 2h = A_1+A_2+A_3 \\ 0 = -\dfrac h2 A_1+\dfrac h2 A_3 \\ \dfrac{2h^3}{3} = \dfrac{h^2}{4}A_1+\dfrac{h^2}{4}A_3 \end{cases}

解得：

A1=A3=4h3,A2=−2h3A_1=A_3=\frac{4h}{3},\qquad A_2=-\frac{2h}{3}

从而求积公式为：

∫−hhf(x) dx≈4h3f ⁣(−h2)−2h3f(0)+4h3f ⁣(h2)\int_{-h}^{h}f(x)\,dx \approx\frac{4h}{3}f\!\left(-\frac{h}{2}\right)-\frac{2h}{3}f(0)+\frac{4h}{3}f\!\left(\frac{h}{2}\right)

对于 f(x)=x3f(x)=x^3， 左端 =0=0， 右端 =4h3(−h38)+4h3(h38)=0=\dfrac{4h}{3}\left(-\dfrac{h^3}{8}\right)+\dfrac{4h}{3}\left(\dfrac{h^3}{8}\right)=0 故成立！

对于 f(x)=x4f(x)=x^4， 左端 =2h55=\dfrac{2h^5}{5}， 右端 =4h3⋅h416⋅2=h56≠2h55=\dfrac{4h}{3}\cdot\dfrac{h^4}{16}\cdot 2=\dfrac{h^5}{6}\neq\dfrac{2h^5}{5} 不成立。

从而代数精度为 33。

根据课上定理，余项形式为：

R2(f)=Ch5f(4)(ξ),ξ∈(−h, h)R_2(f) = C h^5 f^{(4)}(\xi),\qquad \xi\in(-h,\,h)

取 f(x)=x4f(x)=x^4 确定误差常数：

R2(x4)=2h55−h56=7h530R_2(x^4) =\frac{2h^5}{5}-\frac{h^5}{6} =\frac{7h^5}{30}

又因为 f(4)(x)≡24f^{(4)}(x) \equiv 24，可得 C=7720C=\dfrac{7}{720}。

因此截断误差为：

R2(f)=7h5720 f(4)(ξ),ξ∈(−h, h)R_2(f)=\frac{7h^5}{720}\,f^{(4)}(\xi),\qquad \xi\in(-h,\,h)

题 7.4

注意到 Lagrange 插值基函数都是次数不超过 n−1n-1 的多项式，且求积公式代数精确度不小于 n−1n-1。 从而代入 f(x)=li(x),  i=1,⋯ ,nf(x)=l_i(x),\; i=1, \cdots, n 得：

∫abli(x) dx=∑k=1nAkli(xk)=Ak\int_a^b l_i(x)\,dx = \sum_{k=1}^{n}A_k l_i(x_k) = A_k

也即

Ak=∫ablk(x) dx,k=1,2,…,nA_k = \int_a^b l_k(x)\,dx ,\qquad k=1,2,\dots,n

题 7.5

给定 1<p≤100001<p\le 10000，需计算 pp-范数下的圆周率：

πp=2p∫01[u1−p+(1−u)1−p]1/pdu\pi_p=\frac{2}{p}\int_0^1\Bigl[u^{1-p}+(1-u)^{1-p}\Bigr]^{1/p}\mathrm{d}u

记被积函数：

f(u)=[u1−p+(1−u)1−p]1/pf(u)=\bigl[u^{1-p}+(1-u)^{1-p}\bigr]^{1/p}

当 u→0+u\to 0^+ 时，(1−u)1−p→1(1-u)^{1-p}\to 1，故：

f(u)∼u(1−p)/p=u−μ,μ=p−1p∈(0,1)f(u)\sim u^{(1-p)/p}=u^{-\mu},\qquad \mu=\frac{p-1}{p}\in(0,1)

由对称性 f(u)=f(1−u)f(u)=f(1-u)，在 u→1−u\to 1^- 处具有相同的奇异性。

转化被积函数

将 f(u)f(u) 的奇异主项提取出来，令 μ=(p−1)/p\mu=(p-1)/p，将 f(u)f(u) 分解为：

f(u)=u−μ+(1−u)−μ+h(u)f(u) = u^{-\mu}+(1-u)^{-\mu}+h(u)

其中 h(u)=f(u)−u−μ−(1−u)−μh(u)=f(u)-u^{-\mu}-(1-u)^{-\mu}。

在 u→0+u\to 0^+ 处，将 f(u)f(u) 提取 u−μu^{-\mu} 后作泰勒展开：

f(u)=u−μ[1+up−1(1−u)1−p]1/p=u−μ+1pu(p−1)2/p(1−u)1−p+⋯f(u) =u^{-\mu}\Bigl[1+u^{p-1}(1-u)^{1-p}\Bigr]^{1/p} =u^{-\mu}+\frac{1}{p}u^{(p-1)^2/p}(1-u)^{1-p}+\cdots

由于 (p−1)2/p>0(p-1)^2/p>0，第二项及更高阶项在 u→0+u\to 0^+ 时均趋于 00。因此：

lim⁡u→0+h(u)=0−1=−1\lim_{u\to 0^+}h(u)=0-1=-1

由对称性同理可得 lim⁡u→1−h(u)=−1\lim_{u\to 1^-}h(u)=-1。 故 h(u)h(u) 在 [0,1][0,1] 上连续有界，无奇异性。

进一步提高光滑性

为进一步提高数值积分效率，利用对称性 h(u)=h(1−u)h(u)=h(1-u)，将积分压缩至左半区间：

∫01h(u) du=2∫01/2h(u) du.\int_0^1 h(u)\,\mathrm{d}u=2\int_0^{1/2} h(u)\,\mathrm{d}u.

再令 u=t2u=t^2，当 u∈[0,1/2]u\in[0,1/2] 时 t∈[0,1/2]t\in[0,1/\sqrt{2}]，得

∫01h(u) du=2∫01/2k(t) dt\int_0^1 h(u)\,\mathrm{d}u=2\int_0^{1/\sqrt{2}} k(t)\,\mathrm{d}t

其中 k(t)=2t⋅h(t2)k(t)=2t\cdot h(t^2)。

积分公式

奇异主项的积分可解析求出：

∫01u−μ du=∫01(1−u)−μ du=11−μ=p\int_0^1 u^{-\mu}\,\mathrm{d}u=\int_0^1(1-u)^{-\mu}\,\mathrm{d}u=\frac{1}{1-\mu}=p

于是原积分化为：

πp=4p∫01/2k(t) dt+4\pi_p=\frac{4}{p}\int_0^{1/\sqrt{2}} k(t)\,\mathrm{d}t + 4

数值积分方法

对 k(t)k(t) 在 [0,1/2][0,1/\sqrt{2}] 上分段使用 5 点 Gauss-Legendre 公式。 因为其代数精度高，收敛速度快，可以充分利用 k(t)k(t) 的光滑性。

同时为了达到合适的精度，从 N=8N=8 个子区间开始，每次加倍，直到相邻两次积分结果之差的绝对值 <10−14<10^{-14}，可稳定保证 10 位有效数字。

题 7.6

优化问题定义

取 NN 个控制点 ri=(xi,yi)\mathbf r_i=(x_i,y_i)，i=0,…,N−1i=0,\dots,N-1，要求 rN=r0\mathbf r_N=\mathbf r_0。 参数区间 [0,2π][0,2\pi] 均匀划分：θi=2πi/N\theta_i=2\pi i/N。

则投影面积近似为：

A(r)=12∑i=0N−1(xiyi+1−xi+1yi)A(\mathbf r)=\frac12\sum_{i=0}^{N-1}(x_i y_{i+1}-x_{i+1}y_i)

记第 ii 段空间向量 si=(Δxi,Δyi,Δhi)\mathbf s_i=(\Delta x_i,\Delta y_i,\Delta h_i)， 其模 ∥si∥=Δxi2+Δyi2+Δhi2\|\mathbf s_i\|=\sqrt{\Delta x_i^2+\Delta y_i^2+\Delta h_i^2}， 其中 Δxi=xi+1−xi\Delta x_i=x_{i+1}-x_i，Δyi=yi+1−yi\Delta y_i=y_{i+1}-y_i，Δhi=h(xi+1,yi+1)−h(xi,yi)\Delta h_i=h(x_{i+1},y_{i+1})-h(x_i,y_i)。

则空间弧长近似为：

L(r)=∑i=0N−1∥si∥L(\mathbf r)=\sum_{i=0}^{N-1} \|\mathbf s_i\|

优化目标即为：

max⁡r∈R2N  A(r)s.t.L(r)=L0\max_{\mathbf r\in\mathbb R^{2N}}\;A(\mathbf r)\quad\text{s.t.}\quad L(\mathbf r)=L_0

采用罚函数法将约束嵌入优化目标：

max⁡r  Φμ(r)=A(r)−μ2(L(r)−L0)2\max_{\mathbf r}\;\Phi_\mu(\mathbf r)=A(\mathbf r)-\frac{\mu}{2}\bigl(L(\mathbf r)-L_0\bigr)^2

μ>0\mu>0 为罚因子，迭代中逐步增大（μk+1=cμk,  c>1\mu_{k+1}=c\mu_k,\;c>1），驱使约束趋于精确满足。

梯度与步长计算

面积梯度：

∂A∂xi=12(yi+1−yi−1),∂A∂yi=12(xi−1−xi+1)\frac{\partial A}{\partial x_i}=\frac12(y_{i+1}-y_{i-1}),\qquad \frac{\partial A}{\partial y_i}=\frac12(x_{i-1}-x_{i+1})

其中下标模 NN 运算，即 x−1=xN−1x_{-1}=x_{N-1}，xN=x0x_N=x_0。

弧长梯度：

由链式法则：

∂L∂xi=−Δxi+Δhi⋅hx(ri)∥si∥+Δxi−1+Δhi−1⋅hx(ri)∥si−1∥\frac{\partial L}{\partial x_i} =-\frac{\Delta x_i+\Delta h_i\cdot h_x(\mathbf r_i)}{\|\mathbf s_i\|} +\frac{\Delta x_{i-1}+\Delta h_{i-1}\cdot h_x(\mathbf r_i)}{\|\mathbf s_{i-1}\|} ∂L∂yi=−Δyi+Δhi⋅hy(ri)∥si∥+Δyi−1+Δhi−1⋅hy(ri)∥si−1∥\frac{\partial L}{\partial y_i} =-\frac{\Delta y_i+\Delta h_i\cdot h_y(\mathbf r_i)}{\|\mathbf s_i\|} +\frac{\Delta y_{i-1}+\Delta h_{i-1}\cdot h_y(\mathbf r_i)}{\|\mathbf s_{i-1}\|}

其中 hx,hyh_x, h_y 采用中心差商计算：

hx(xi,yi)≈h(xi+δ,yi)−h(xi−δ,yi)2δhy(xi,yi)≈h(xi,yi+δ)−h(xi,yi−δ)2δ\begin{aligned} h_x(x_i,y_i)&\approx\frac{h(x_i+\delta,y_i)-h(x_i-\delta,y_i)}{2\delta} \\ h_y(x_i,y_i)&\approx\frac{h(x_i,y_i+\delta)-h(x_i,y_i-\delta)}{2\delta} \end{aligned}

目标函数梯度：

∇Φμ=∇A−μ(L−L0)∇L\nabla\Phi_\mu=\nabla A-\mu(L-L_0)\nabla L

线搜索确定步长：

从步长 α=α0\alpha=\alpha_0 开始，每次折半，直到

Φμ(r+αg)≥Φμ(r)+σα gTg\Phi_\mu(\mathbf r+\alpha\mathbf g)\geq\Phi_\mu(\mathbf r)+\sigma\alpha\,\mathbf g^{\mathsf T}\mathbf g

其中取 σ=0.1\sigma = 0.1。也即取步长使得提升达到线性情况的 σ\sigma 倍。

整体算法

算法步骤：

  1. 初始化：给定闭合控制点 r(0)\mathbf{r}^{(0)}，设置罚因子 μ=μ0\mu=\mu_0。
  2. 外循环：若 ∣L(r)−L0∣<δ|L(\mathbf{r})-L_0|<\delta 则结束；否则进入内循环。
  3. 内循环：
    • 计算所有点的高度 h(ri)h(\mathbf{r}_i) 及地形梯度 hx,hyh_x,h_y；
    • 计算面积梯度 ∇A\nabla A 与弧长梯度 ∇L\nabla L；
    • 组装 g=∇Φμ=∇A−μ(L−L0)∇L\mathbf{g}=\nabla\Phi_\mu=\nabla A-\mu(L-L_0)\nabla L；
    • 线搜索确定步长 α\alpha；
    • 更新 r←r+αg\mathbf{r}\leftarrow\mathbf{r}+\alpha\mathbf{g}，其中强制 rN=r0\mathbf{r}_N=\mathbf{r}_0；
    • 若 ∥g∥<ε\|\mathbf{g}\|<\varepsilon 则退出内循环，否则重复。
  4. 增大罚因子：μ←cμ\mu\leftarrow c\mu，返回步骤 2。
  5. 输出：r∗\mathbf{r}^* 与最大面积 A∗A^*。

复杂度分析

操作复杂度
面积 AA 及其梯度O(N)O(N)
弧长 LL 及其梯度O(N)O(N)
每次梯度上升迭代O(N)O(N)

总体复杂度分析取决于内外层迭代次数。