PTA | 数值分析 | 题解汇总

【6-1】Numerical Summation of a Series

求：$\varphi (x)$，其中 $x$ 取 $0.0,0.1,0.2,...,300.0$
条件：
1）$\varphi (x)={\sum }_{k=1}^{\infty }\frac{1}{k(k+x)}$
2）对于输出的每一项，要求绝对误差小于 $10^{-10}$
3）时限为 $100$ ms

分析：一种简单暴力的想法是设置一个阈值 $T$：对 $k\le T$，直接暴力累加；对 $k>T$，直接估计。然而，在指定的时间内，不管如何设置 $T$，精度都远远不够。

考虑对 $\varphi (x)$ 做一些转化，使得转化后的对象 $\mu (x)$ 收敛得更快，对转化后的对象 $\mu (x)$ 用相同的方法计算，然后再根据 $\varphi (x)$ 与 $\mu (x)$ 的关系，计算得到 $\varphi (x)$。

注意到
$$\varphi (1)=\sum _{k=1}^{\infty }\frac{1}{k(k+1)}=1$$
对 $\varphi (x)$ 作转化：令 $\varphi (x)$ 减去 $\varphi (1)$，得
$$\varphi (x)-\varphi (1)=(1-x)\sum _{k=1}^{\infty }\frac{1}{k(k+1)(k+x)}=(1-x)\mu (x)$$
然而这样还是过不了。

正解的做法应该是对 $\mu (x)$ 继续进行变换，然而我又变换了三四次，依然无法满足题意要求，误差还越来越大（我猜测这是因为对转化量的的要求越来越高）。

只能去考虑对于不同的 $x$，$\varphi (x)$ 之间存在什么联系。

先对 $\varphi (x)$ 进行裂项转化：
$$\varphi (x)=\sum _{k=1}^{\infty }\frac{1}{k(k+x)}=\frac{1}{x}\sum _{k=1}^{\infty }(\frac{1}{k}-\frac{1}{k+x})$$

$$\varphi (x+1)=\sum _{k=1}^{\infty }\frac{1}{k(k+x+1)}=\frac{1}{x+1}\sum _{k=1}^{\infty }(\frac{1}{k}-\frac{1}{k+x+1})$$

那么有
$$x\varphi (x)=(x+1)\varphi (x+1)-\frac{1}{1+x}$$
观察得知 $\varphi$ 越小，我们的方法计算精度越高。于是我们可以设置更大的阈值 $T$，更精确地计算 $\varphi (0.0,0.1,0.2,...,0.9)$，然后根据上面的等式关系，推出所有的解。

---

【6-2】Root of a Polynomial

求：$x$
已知：
1）$n$ 次多项式 $p(x)=c_n{x}^n+c_{n-1}{x}^{n-1}+...+c_{1}x+c_0$
2）$a,b$
条件：
1）$p(x)=0$
2）$a<x<b$
2）保证 $x$ 有且仅有一个解

注意事项：
1）有一个测试点在 $x=b$ 处有零点，但实际答案为 $(a,b)$ 中的另一个零点。
2）在做牛顿迭代法的时候，需要用一个更优秀的迭代公式，以及需要取多个迭代的初始值。

---

【6-3】There is No Free Lunch
求：$c_i$
已知：$n(\le 10000)$, $p_1,p_2,...,p_n$
条件：
$$\left(\begin{array}{cccccccc}2& \frac{1}{2}& & & & & & \frac{1}{2}\\ \frac{1}{2}& 2& \frac{1}{2}& & & & & \\ & \frac{1}{2}& 2& \frac{1}{2}& & & & \\ & & \frac{1}{2}& 2& \frac{1}{2}& & & \\ & & & \ddots & \ddots & \ddots & & \\ & & & & \frac{1}{2}& 2& \frac{1}{2}& \\ & & & & & \frac{1}{2}& 2& \frac{1}{2}\\ \frac{1}{2}& & & & & & \frac{1}{2}& 2\end{array}\right)\left(\begin{array}{c}{c}_1\\ c_2\\ c_3\\ c_4\\ ⋮\\ c_{n-2}\\ c_{n-1}\\ c_n\end{array}\right)=\left(\begin{array}{c}{p}_1\\ p_2\\ p_3\\ p_4\\ ⋮\\ p_{n-2}\\ p_{n-1}\\ p_n\end{array}\right)$$

没有学标准做法，直接自己胡了一个线性算法过了。
用 (1, 1) 消去 (2, 1) 和 (n, 1)。
用 (2, 2) 消去 (3, 2) 和 (n, 2)。
用 (3, 3) 消去 (4, 3) 和 (n, 3)。
…
用 (n - 2, n - 2) 消去 (n - 1, n - 2) 和 (n - 2, n - 2)。
可以发现，我们只需要维护主对角线的元素、最后一列的元素以及该行元素线性组合得到的值，并得到如下矩阵：
$$\left(\begin{array}{cccccccc}?& ?& & & ...& & ?& B_1\\ & ?& ?& & ...& & ?& B_2\\ & & ?& ?& ...& & ?& B_3\\ & & & \ddots & \ddots & & ⋮& ⋮\\ & & & & ?& ?& ?& B_{n-2}\\ & & & & & ?& ?& B_{n-1}\\ & & & & & ?& ?& B_n\end{array}\right)$$
首先，我们用行列式解出 $c_{n-1},c_n$，然后可以轻松递推得到 $c_{n-2},c_{n-3},...,c_1$。

---

【6-4】Compare Methods of Jacobi with Gauss-Seidel

先调整主元，然后就是模板题。

Jacobi：

for (int it = 1; it <= MAXN; it ++) {
	for (int i = 0; i < n; i ++) {
		_x[i] = b[i];
		for (int j = 0; j < n; j ++)
			if (i != j)
				_x[i] -= a[i][j] * x[j];
		_x[i] /= a[i][i];
	}
}

Gauss：

for (int it = 1; it <= MAXN; it ++) {
	for (int i = 0; i < n; i ++) {
		_x[i] = b[i];
		for (int j = 0; j < i; j ++)
			_x[i] -= a[i][j] * _x[j];
	for (int j = i + 1; j < n; j ++)
		_x[i] -= a[i][j] * x[j];
		_x[i] /= a[i][i];
}

---

【6-5】Approximating Eigenvalues

对于 $n$ 阶矩阵 $A$，如何求 $A$ 的最大特征值 ${\lambda }_1$？
直接有
$${\lambda }_1\approx \frac{(A^k{)}_{i,j}}{(A^{k-1}{)}_{i,j}}$$

对于 $n$ 阶矩阵 $A$，如何求 $A$ 的最小特征值 ${\lambda }_n$？
先求 $A^{-1}$ 的最大特征值 ${\lambda }_1^{\prime }$，则有 ${\lambda }_n=1/{\lambda }_1^{\prime }$。

对于 $n$ 阶矩阵 $A$，给定某个近似特征值 ${\lambda }_0$ 以及近似特征向量 $x_0$，如何更加精确地逼近？
$\lambda -{\lambda }_0$ 是矩阵 $A-{\lambda }_{0}I$ 的最小特征值。

在解最小特征值的时候，我的做法是直接将矩阵的逆 $A^{-1}$ 求出来，然后每次暴力乘上 $A^{-1}$ 算特征值，最后取反。

---

【6-6】Cubic Spline

直接高斯消元法暴力解方程可过。为了建立方程组的简洁，这是我定义的宏：

#define M (4 * N)
#define A(x) (4 * (x) - 3)
#define B(x) (4 * (x) - 2)
#define C(x) (4 * (x) - 1)
#define D(x) (4 * (x) - 0)
#define NEW (tot ++)
#define ADD(x, y) (eq[tot][(x)] = (y))
#define LIN(i) (x[i] - x[i - 1])
#define SQR(i) (LIN(i) * LIN(i))
#define CUB(i) (SQR(i) * LIN(i))

这样就可以把建立方程组写得很简洁：

if (Type == 1) {
		NEW, ADD(B(1), 1), ADD(M + 1, S0);
		NEW, ADD(B(N), 1), ADD(C(N), 2 * LIN(N)), ADD(D(N), 3 * SQR(N)), ADD(M + 1, SN);
	}
	else {
		NEW, ADD(C(1), 2), ADD(M + 1, S0);
		NEW, ADD(C(N), 2), ADD(D(N), 6 * LIN(N)), ADD(M + 1, SN);
	}
	for (int I = 0; I <= N; I ++) {
		if (I > 0)
			NEW, ADD(A(I), 1), ADD(B(I), LIN(I)), ADD(C(I), SQR(I)), ADD(D(I), CUB(I)), ADD(M + 1, f[I]);
		if (I < N)
			NEW, ADD(A(I + 1), 1), ADD(M + 1, f[I]);
	}
	for (int I = 1; I <= N - 1; I ++) {
		NEW, ADD(B(I), 1), ADD(C(I), 2 * LIN(I)), ADD(D(I), 3 * SQR(I)), ADD(B(I + 1), -1);
		NEW, ADD(C(I), 2), ADD(D(I), 6 * LIN(I)), ADD(C(I + 1), -2);
	}

---

【6-7】Orthogonal Polynomials Approximation

在实现上，我不想写类，直接用数组存储多项式。
写得有点麻烦，可以感受一下：

	int k = 1;
	while (k < MAX_n && fabs(err) >= *eps) {
		k ++;
		cpy(tmp, phi[k - 1]);
		mov(tmp);
		double bk = inner2(tmp, phi[k - 1]) / inner2(phi[k - 1], phi[k - 1]);
		double ck = inner2(tmp, phi[k - 2]) / inner2(phi[k - 2], phi[k - 2]);
		cpy(phi[k], phi[k - 1]);
		mov(phi[k]);
		cpy(tmp, phi[k - 1]);
		mul(tmp, -bk);
		add(phi[k], tmp);
		cpy(tmp, phi[k - 2]);
		mul(tmp, -ck);
		add(phi[k], tmp);
		a[k] = inner1(phi[k], f) / inner2(phi[k], phi[k]);
		cpy(tmp, phi[k]);
		mul(tmp, a[k]);
		add(c, tmp);
		err -= a[k] * inner1(phi[k], f);
	}

---

【6-8】Shape Roof

微元法分析：

$$\begin{gathered} dl=dx\ast \sqrt{1+{\frac{dy}{dx}}^2} \\ l={\int }_a^b\sqrt{1+{\frac{dy}{dx}}^2}dx={\int }_a^b\text{f0}(x)dx \end{gathered}$$

然后直接模板题。
注意设 EPS 为 1e-9，比 1e-9 大则精度不够答案错误，比 1e-9 小则不能计算出来死循环。