给我一个使用基于Galerkin原理的arnoldi算法迭代求解线性方程组的例子

考虑如下线性方程组：

$$Ax=b$$

其中，$A$是一个$n\times n$的稀疏矩阵，$b$是一个$n\times 1$的向量，我们需要求解$x$。

使用Galerkin原理，我们可以将该线性方程组转化为求解下面的矩阵特征值问题：

$$H y_k = y_{k+1}$$

其中，$H=V^T A V$，$V$是一个$n\times k$的正交矩阵，$y_k$是一个$k\times 1$的向量，$y_{k+1}$是$y_k$在$V$的列空间上的投影。

Arnoldi算法是一种迭代方法，用于求解上面的矩阵特征值问题。它的基本思想是在每一步迭代中，将当前的正交基向量$v_k$与$A$相乘，得到一个新的向量$h_{k+1,k}=Av_k$，然后对$h_{k+1,k}$进行正交化，得到一个新的正交基向量$v_{k+1}$，以此类推。

具体来说，Arnoldi算法的迭代过程如下：

1. 选择一个初始向量$v_1$，将其单位化得到$v_1^*$。

2. 对于$k=1,2,\cdots,m$，做如下迭代：

(a) 计算$h_{k+1,k}=Av_k$。

(b) 令$H_{k+1,k}=h_{1:k+1,k}$，$V_{k+1}=[v_1,v_2,\cdots,v_{k+1}]$。

(c) 对$H_{k+1,k}$进行正交化，得到一个新的向量$h_{k+1,k}^*$，并令$v_{k+1}=h_{k+1,k}^*- \sum_{i=1}^k h_{i,k}^*v_i$。

(d) 将$v_{k+1}$单位化得到$v_{k+1}^*$。

3. 将$V_m$和$H_{m,m-1}$代入矩阵特征值问题$H y_k = y_{k+1}$中，使用标准的特征值求解方法求解特征值和特征向量。

4. 对于求得的特征值$\lambda_k$和特征向量$y_k$，设$x_k=V_my_k$，则$x_k$是原线性方程组的一个解。

下面给出一个具体的例子：

考虑如下线性方程组：

$$\begin{pmatrix}1 & 2 & 0 \\ 2 & 3 & 1 \\ 0 & 1 & 4\end{pmatrix}\begin{pmatrix}x_1 \\ x_2 \\ x_3\end{pmatrix}=\begin{pmatrix}1 \\ 2 \\ 3\end{pmatrix}$$

我们使用Arnoldi算法迭代求解该线性方程组。首先选取一个单位向量$v_1=[1,0,0]^T$，然后进行迭代。迭代过程中，我们需要计算$Av_k$，并对其进行正交化，得到一个新的向量$v_{k+1}$。

迭代1：

$$h_{2,1}=Av_1=\begin{pmatrix}1 \\ 2 \\ 0\end{pmatrix}$$

进行正交化，得到$h_{2,1}^*=[\frac{1}{\sqrt{5}},\frac{2}{\sqrt{5}},0]^T$，然后令$v_2=h_{2,1}^*$，得到$v_2^*=[\frac{1}{\sqrt{5}},\frac{2}{\sqrt{5}},0]^T$。

迭代2：

$$h_{3,2}=Av_2=\begin{pmatrix}\frac{2}{\sqrt{5}} \\ \frac{7}{\sqrt{5}} \\ \frac{1}{\sqrt{5}}\end{pmatrix}$$

进行正交化，得到$h_{3,2}^*=[\frac{2}{\sqrt{30}},\frac{7-2\sqrt{5}}{\sqrt{150}},\frac{1}{\sqrt{30}}]^T$，然后令$v_3=h_{3,2}^*-\frac{h_{2,2}^*}{\|h_{2,2}^*\|^2}(h_{3,2}^*\cdot h_{2,2}^*)$，其中，$h_{2,2}=v_2^* A v_2=\begin{pmatrix}\frac{2}{\sqrt{5}} & \frac{7}{\sqrt{5}} & 0\end{pmatrix}^T$。化简得：

$$v_3^*=[\frac{2}{\sqrt{30}},\frac{3\sqrt{5}-7}{\sqrt{150}},\frac{1}{\sqrt{30}}]^T$$

迭代3：

此时，$H_{4,3}$已经是一个上Hessenberg矩阵，无需进行正交化。我们可以将$V_3$和$H_{3,2}$代入矩阵特征值问题$H y_k = y_{k+1}$中，使用标准的特征值求解方法求解特征值和特征向量。求解得到的特征值为$\lambda_1\approx 4.5804$，$\lambda_2\approx 1.7484$，$\lambda_3\approx 0.6712$，对应的特征向量为：

$$y_1^*=[0.1708,-0.9849,-0.0224], y_2^*=[0.3390,-0.1434,0.9294], y_3^*=[-0.9246,-0.0895,0.3702]$$

最终，我们得到的三个解分别为$x_1=V_3y_1=[-0.0150,-1.2553,0.3243]^T$，$x_2=V_3y_2=[-0.9265,0.8456,0.9272]^T$，$x_3=V_3y_3=[-0.6292,0.3155,1.2036]^T$。