科学计算作业 三

科学计算

陆一洲 5140309557

3.30

考察要选取的元素的符号，以避免它被消。
若为正，则 $\alpha$ 为负，反之则为正。

3.33

(a) $LU$ 矩阵可以同时存在 $A$ 中，需要 $1\times n$ 的空间存储置换。
(b) 需要 $1$ 个 $m\times m$ 矩阵分别记录当前 $Q$ 和 $1$ 个 $1\times m$ 的向量存储 $v$。矩阵 $H$ 由 $v$ 隐式地存储，消元后的 $A$ 即 $R$。

3.44

法一：高斯消元。在选取列主元时，若所有元素的绝对值都很小，则记该列元素为 $0$，矩阵中非零行的个数即为矩阵的秩。
法二：计算奇异值矩阵，绝对值大于某个极小值的元素个数即为矩阵的秩。

3.45

(a) $$mn-\frac{1}{3}n^3=\frac{5}{3}n^3$$
(b) $$\frac{1}{2}mn-\frac{1}{6}n^3=\frac{11}{6}n^3$$
(c) 速度排序：norm equation > QR 分解 > 奇异值分解。

3.1

(a) $$\begin{pmatrix} 1&10 \\ 1&15 \\ 1&20 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}= \begin{pmatrix} 11.60 \\ 11.85 \\ 12.25 \end{pmatrix}$$
(b) 无解
(c)
$$A^TA= \begin{pmatrix} 3&45 \\ 45&725 \end{pmatrix} ,\ \ \ A^Tb= \begin{pmatrix} 35.7 \\ 538.75 \end{pmatrix} \\ A^TAx=A^Tb\ \ \ \Rightarrow\ \ \ x_1=10.925,2_2=0.065 \\ Ax= \begin{pmatrix} 11.575\\ 11.9\\ 12.225 \end{pmatrix} , \ \ \ \Delta b= \begin{pmatrix} -0.025 \\ 0.05\\ -0.025 \end{pmatrix}$$

3.3

$$\begin{pmatrix} 1 & 2.71828\\ 2 & 7.38906\\ 3 & 20.0855 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} =\begin{pmatrix} 2\\3\\5 \end{pmatrix}\\ \begin{pmatrix} 14 &77.7529 \\ 77.7529 &465.415 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} =\begin{pmatrix} 23 \\ 128.031 \end{pmatrix}\\ x_1=1.59423,\ x_2=0.00875469$$ 验证：
$$Ax=\begin{pmatrix} 1.61803 \\ 3.24315\\ 4.95753 \end{pmatrix}$$ 与 $b$ 大致相同。

3.14

$$\begin{align*} H^2&=(I-2\frac{vv^T}{v^Tv})^2\\ &=I-4\frac{vv^T}{v^Tv}+4\frac{(vv^T)^2}{v^T(vv^T)v}\\ &=I-4\frac{vv^T}{v^Tv}+4\frac{vv^T}{v^Tv}\\ &=I\\ H^T&=(I-2\frac{vv^T}{v^Tv})^T\\ &=I-(2\frac{vv^T}{v^Tv})^T\\ &=I-2\frac{vv^T}{v^Tv}\\ &=H \end{align*}\\ $$ 所以 $H$ 是正交且对称的。

3.31

下证：当 $x=\sum_{\sigma_i\neq 0}\frac{u_i^Tb}{\sigma_i}v_i$ 时，有 $$A^TAx=A^Tb$$ 证明： $$\begin{align*} A^TAx&=A^TA\sum_{\sigma_i\neq 0}\frac{u_i^Tb}{\sigma_i}v_i\\ &=A^TAV\Sigma_1^{-1}U_1^Tb\\ &=V\Sigma^T\Sigma\Sigma_1^{-1}U_1^Tb \end{align*}$$ 其中， $$\Sigma^T\Sigma\Sigma_1^{-1}=\Sigma_1$$ 所以， $$A^TAx=V\Sigma_1U_1^Tb=A^Tb$$ 命题得证。

3.32

(a) $$\begin{align*} AA^+A &= U\Sigma V^TV\Sigma^+U^TU\Sigma V^T\\ &= U\Sigma \Sigma^+\Sigma V^T\\ &= U\Sigma V^T\\ &= A \end{align*}$$ (b) $$\begin{align*} A^+AA^+ &= V\Sigma^+U^TU\Sigma V^TV\Sigma^+U^T\\ &= V\Sigma^+\Sigma \Sigma^+U^T\\ &= V\Sigma^+U^T\\ &= A^+ \end{align*}$$ (c) $$\begin{align*} (AA^+)^T &= (U\Sigma V^TV\Sigma^+U^T)^T \\ &= U(\Sigma\Sigma^+)^TU^T \\ &= U\Sigma\Sigma^+U^T \\ &= U\Sigma V^TV\Sigma^+U^T \\ &= AA^+ \end{align*}$$ (d) $$\begin{align*} (A^+A)^T &= (V\Sigma^+U^TU\Sigma V^T)^T \\ &= V(\Sigma^+\Sigma)^TV^T \\ &= V\Sigma^+\Sigma V^T \\ &= V\Sigma^+U^TU\Sigma V^T \\ &= A^+A \end{align*}$$

编程题

核心代码如下：

Matrix A(m,n);
vector<double> b(m), idx(n);
for (int i=0; i<m; i++){
    double tmp=Git<double>();
    b[i]=Git<double>();
    A[i][0]=1;
    for (int j=1; j<n; j++)
        A[i][j]=A[i][j-1]*tmp;
}
for (int i=0; i<n; i++) idx[i]=i;
Matrix A0(A);
vector<double> b0(b), x(M);
Householder_Transformation(FLAG);
for (int i=n-1; i>=0; i--){
    double tmp=b[i];
    for (int j=i+1; j<n; j++)
        tmp-=A[i][j]*x[j];
    x[i]=t/A[i][i];
}
vector<double> x0(x);
for (int i=0; i<n; i++) x[idx[i]]=x0[i];
for (int i=0; i<m; i++){
    b[i]=0;
    for (int j=0; j<n; j++)
        b[i]+=A0[i][j]*x[j];
        b[i]=b0[i]-b[i];
}
cout<<Norm(b)<<endl;

记 $\tau=\|b-Ax\|_2$。

• 没有主元操作时， $\bar\tau=1.45e-7$；
• 有主元操作时，$\bar\tau=9.93e-6$；

有主元后误差更大，其原因为对于 $x_i,x_j\ (i<j)$，在 $A$ 中 $x_i$ 的系数比 $x_j$ 大若干倍。所以计算时 $x_i$ 产生的误差更大；
有主元操作时，可以先计算 $x_j$ 再计算 $x_i$。此时，$x_i$ 的误差会由于之前计算的累积误差被进一步放大，导致拟合结果变得更差。