【6-1】Numerical Summation of a Series
求：\(\phi(x)\)，其中 \(x\) 取 \(0.0, 0.1, 0.2, ..., 300.0\)
条件：
1）\(\phi(x) = \sum_{k = 1} ^ {\infty} \frac{1}{k(k + x)}\)
2）对于输出的每一项，要求绝对误差小于 \(10 ^ {-10}\)
3）时限为 \(100\) ms
分析：一种简单暴力的想法是设置一个阈值 \(T\)：对 \(k \le T\)，直接暴力累加；对 \(k > T\)，直接估计。然而，在指定的时间内，不管如何设置 \(T\)，精度都远远不够。
考虑对 \(\phi(x)\) 做一些转化，使得转化后的对象 \(\mu(x)\) 收敛得更快，对转化后的对象 \(\mu(x)\) 用相同的方法计算，然后再根据 \(\phi(x)\) 与 \(\mu(x)\) 的关系，计算得到 \(\phi(x)\)。
注意到
对 \(\phi(x)\) 作转化：令 \(\phi(x)\) 减去 \(\phi(1)\)，得
然而这样还是过不了。
正解的做法应该是对 \(\mu(x)\) 继续进行变换，然而我又变换了三四次，依然无法满足题意要求，误差还越来越大（我猜测这是因为对转化量的的要求越来越高）。
只能去考虑对于不同的 \(x\)，\(\phi(x)\) 之间存在什么联系。
先对 \(\phi(x)\) 进行裂项转化：
那么有
观察得知 \(\phi\) 越小，我们的方法计算精度越高。于是我们可以设置更大的阈值 \(T\)，更精确地计算 \(\phi(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_1x + 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\)
条件：
没有学标准做法，直接自己胡了一个线性算法过了。
用 (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)。
可以发现，我们只需要维护主对角线的元素、最后一列的元素以及该行元素线性组合得到的值，并得到如下矩阵：
首先，我们用行列式解出 \(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\)？
直接有
对于 \(n\) 阶矩阵 \(A\)，如何求 \(A\) 的最小特征值 \(\lambda_n\)？
先求 \(A ^ {-1}\) 的最大特征值 \(\lambda_1 '\)，则有 \(\lambda_n = 1 / \lambda_1'\)。
对于 \(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
微元法分析：
然后直接模板题。
注意设 EPS 为 1e-9，比 1e-9 大则精度不够答案错误，比 1e-9 小则不能计算出来死循环。