PTA | 数值分析 | 题解汇总

2022-04-01 14:01:29 阅读：192 来源： 互联网

标签：tmp phi frac 题解 over 汇总 PTA ADD define

【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(1) = \sum_{k = 1} ^ {\infty} \frac{1}{k(k + 1)} = 1 \]

对 \(\phi(x)\) 作转化：令 \(\phi(x)\) 减去 \(\phi(1)\)，得

\[\phi(x) - \phi(1) = (1 - x)\sum_{k = 1} ^ {\infty} \frac{1}{k(k + 1)(k + x)} = (1 - x) \mu(x) \]

然而这样还是过不了。

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

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

先对 \(\phi(x)\) 进行裂项转化：

\[\phi(x) = \sum_{k = 1} ^ {\infty} \frac{1}{k(k + x)} = \frac{1}{x} \sum_{k = 1} ^ {\infty} (\frac{1}{k} - \frac{1}{k + x}) \]

\[\phi(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\phi(x) = (x + 1)\phi(x + 1) - \frac{1}{1 + 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\)
条件：

\[\begin{pmatrix} 2 & {1 \over 2} & & & & & & {1 \over 2} \\ {1 \over 2} & 2 & {1 \over 2} & & & \\ & {1 \over 2} & 2 & {1 \over 2} & \\ & & {1 \over 2} & 2 & {1 \over 2} \\ & & & \ddots & \ddots & \ddots \\ & & & & {1 \over 2} & 2 & {1 \over 2} \\ & & & & & {1 \over 2} & 2 & {1 \over 2} \\ {1 \over 2} & & & & & & {1 \over 2} & 2 \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \\ c_3 \\ c_4 \\ \vdots \\ c_{n - 2} \\ c_{n - 1} \\ c_n \end{pmatrix} = \begin{pmatrix} p_1 \\ p_2 \\ p_3 \\ p_4 \\ \vdots \\ p_{n - 2} \\ p_{n - 1} \\ p_n \end{pmatrix} \]

没有学标准做法，直接自己胡了一个线性算法过了。
用 (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)。
可以发现，我们只需要维护主对角线的元素、最后一列的元素以及该行元素线性组合得到的值，并得到如下矩阵：

\[\begin{pmatrix} ? & ? & & & ... & & ? & B_1\\ & ? & ? & & ... & & ? & B_2 \\ & & ? & ? & ... & & ? & B_3 \\ & & & \ddots & \ddots & &\vdots & \vdots \\ & & & & ? & ? & ? & B_{n - 2} \\ & & & & & ? & ? & B_{n - 1} \\ & & & & & ? & ? & B_n \end{pmatrix} \]

首先，我们用行列式解出 \(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 {(A ^ k)_{i, j} \over (A ^ {k - 1})_{i, j}} \]

对于 \(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

微元法分析：

\[dl = dx * \sqrt{1 + {dy \over dx} ^ 2} \\ l = \int_{a} ^ b \sqrt{1 + {dy \over dx} ^ 2} dx = \int_{a} ^ b \text{f0}(x) dx \]

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

标签：tmp,phi,frac,题解,over,汇总,PTA,ADD,define
来源： https://www.cnblogs.com/Sdchr/p/16086609.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

PTA | 数值分析 | 题解汇总