北邮数值与符号计算实验线性代数方程组求解.docx-资料库

fzu031002319-9532742-4744300845177479850.docx.pdf-第1页.png

第1页 / 共6页

fzu031002319-9532742-4744300845177479850.docx.pdf-第2页.png

第2页 / 共6页

fzu031002319-9532742-4744300845177479850.docx.pdf-第3页.png

第3页 / 共6页

fzu031002319-9532742-4744300845177479850.docx.pdf-第4页.png

第4页 / 共6页

fzu031002319-9532742-4744300845177479850.docx.pdf-第5页.png

第5页 / 共6页

fzu031002319-9532742-4744300845177479850.docx.pdf-第6页.png

第6页 / 共6页

1 实验要求用 C/C++语言实现如下函数： bool lu(double* a, int* pivot, int n); 矩阵的 LU 分解 bool guass(double const* lu, int const* p, double* b, int n); 求线性代数方程组的解 void qr(double* a, double* d, int n); 矩阵的 QR 分解 bool hshld(double const*qr, double const*d, double*b, int n); 求线性代数方程组的解用上面的函数求解下面的问题： 1.1 求解 n 阶 Hilbert 方阵的线代数方程组 Hilbert 矩阵在第 j 行，第 k 列处的值为 = 1 +−1 适当设置 b 的值使得理论解围 x=(1,1,…,1)T。列选主元 Gauss 消去法与 Householder 变换法分别对 n=10,20,30 求解。给出数值解与理论解之差的二范数。 1.2 求 Gn 矩阵的线代数方程组用两种方法对 n=50,100,150 求解。理论值为 x=(1,1,…,1)T。给出数值解与理论解之差的二范数。令 n×n 矩阵 0 0 1 令 n 维向量= 2,1,0,−1,−2,…,5−,4−,2− 再令= 2,3,4…,, 用两种方法对 n=50,100,150 求解。给出数值解的残量的二范数。 0 ⋯ 0 0 1 0 ⋯ 0 1 −1 0 ⋯ 0 −1 −1 ⋯ −1 −1 −1 −1 ⋯ 1 −1 −1 −1 −1 ⋯ −1 = 1 1 1 1 1 2 实验平台 IDE：Visual Studio 2013 语言：C++ 操作系统：Windows 10 3 函数实现 3.1 矩阵的 LU 分解 bool lu(double* a, int* pivot, int n){ int i, j, k; double max = 0, temp = 0.0; for (i = 0; i < n - 1; i++){ //对第i列进行处理 // 选出第i列的主元,记录主元位置 max = fabs(a[n*i + i]); pivot[i]=i; for (j = i + 1; j < n; j++){ //第j行

if (fabs(a[n*j + i]) > max) { max = fabs(a[n*j + i]); pivot[i] = j; } } //若最大的列主元为0，则不能进行LU分解 if (fabs(a[n*n - 1]) < 1e-14) { printf("矩阵不能进行LU分解！\n"); return true; } // 对第i列进行行变换，使主元在对角线上 if(pivot[i]!=i) { for (j = i; j < n; j++){ //i,j与pivot[i],j交换对上三角处理 temp = a[n*i + j]; a[n*i + j] = a[n*pivot[i] + j]; a[n*pivot[i] + j] = temp; } } for (j = i + 1; j < n; j++)//Pi 部分下三角L a[n*j + i]=a[n*j + i]/a[n*i + i]; for (j = i + 1; j < n; j++)//计算上三角U for(k=i+1; k

} } //按pivot对b行变换，与LU匹配 for (i = 0; i < n - 1; i++) { temp = b[p[i]]; b[p[i]] = b[i]; b[i] = temp; } //Ly=b,将y的内容放入b for(i=0; i= 0; i--) { for (j = n - 1; j > i; j--) b[i] = b[i] - lu[n*i + j] * b[j]; b[i] = b[i] / lu[n*i + i]; } return false; } 3.3 矩阵的 QR 分解 void qr(double* a, double* d, int n) { int i, j, l, k; double tmp, m; double *temp; temp = (double *)malloc(sizeof(double)*n); for (i = 0; i < n - 1; i++){ //对第i列进行处理 //获得tmp值 m = 0; for (j = i; j < n; j++)//j表示第j行 m = m + a[n*j + i] * a[n*j + i ]; if (a[n*i + i ] > 0) m = -sqrt(m); else m = sqrt(m); //获得temp放入矩阵，并存主元d tmp = 0; d[i] = m; a[n*i + i] = a[n*i + i] - m; for (j = i; j <= n - 1; j++) tmp = tmp + a[n*j + i] * a[n*j + i]; tmp = sqrt(tmp); for (j = i; j <= n - 1; j++) a[n*j + i] = a[n*j + i] / tmp; // 调整矩阵

for (k = i + 1; k < n; k++) { for (j = i; j < n; j++) { tmp = 0; for (l = i; l < n; l++) tmp = tmp + a[n*j + i] * a[n*l + i] * a[n*l + k]; temp[j] = a[j*n + k] - 2 * tmp; } for (j = i; j < n; j++) a[j*n + k] = temp[j]; } } d[n - 1] = a[(n - 1)*n + n - 1]; } 3.4 Householder 变换法求线性代数方程组的解 bool hshld(double const *qr, double const *d, double *b, int n) { int i, j, l; double rem; double *temp; temp = (double *)malloc(sizeof(double)*n); for (i = 0; i < n; i++){ if (fabs(d[i]) < 1e-14) { printf("R矩阵奇异！\n"); return true; } } //计算 b = Qtb for (i = 0; i < n - 1; i++) { for (j = i; j < n; j++) { rem = 0; for (l = i; l < n; l++) rem = rem + qr[l*n + i] * qr[j*n + i] * b[l]; temp[j] = b[j] - 2 * rem; } for (j = i; j < n; j++) b[j] = temp[j]; } //计算 x = R\b for (j = n - 1; j > -1; j--) { for (l = n - 1; l != j; --l) b[j] = b[j] - b[l] * qr[j*n + l]; b[j] = b[j] / d[j]; } return false; }

4 实验结果 4.1 Hilbert 方阵线代数方程组 4.1.1 列选主元 Gauss 消去法 n=10 n=20 n=30 4.1.2Householder 变换法 n=10 n=20 n=30 4.2 n×n 矩阵线代数方程组 4.2.1 列选主元 Gauss 消去法 n=50 n=100 n=150 4.2.2 Householder 变换法 n=50 n=100 n=150 4.3 n×n 矩阵线代数方程组 4.3.1 列选主元 Gauss 消去法 n=50 n=100 n=150 4.3.2 Householder 变换法 n=50 n=100 n=150 5 实验总结从第一个实验的结果可以看到，当 n=10 时，列选主元 Gauss 消去法得到的数值解与理论解之差的二范数，比 Householder 变换法的结果要小，说明此时列选主元 Gauss 消去法更加精确。当 n=20 时，对矩阵进行 LU 分解过程中，最大的列主元绝对值小于 1e-14，认为是 0，已不能进行 LU 分解，所以分解失败，也不能继续求解方程组。Householder 变换法对矩阵进

行 QR 分解后，由于 R 矩阵奇异，所以求解方程组失败。当 n=30 时，结论与 n=20 一致。说明矩阵阶数越大，计算过程中丢失精度也越严重。第二个实验，矩阵元素很简单，主对角元素和最后一列为 1，下三角元素为-1，其余为 0。当 n=50 时，列选主元 Gauss 消去法得到的二范数为 0，而 Householder 变换法得到的二范数虽然很小，但精确度显然不如前者。当 n=100 时，列选主元 Gauss 消去法得到的二范数开始变大，而 Householder 变换法得到的二范数也有所增大，却比前者小得多。当 n=150 时，列选主元 Gauss 消去法得到的二范数继续变大，Householder 变换法得到的二范数也继续增大，依然比前者小得多。说明对于某些矩阵，QR 分解比 LU 更有优势。第二个实验又提出另外一个 bn，当 n=50 时，列选主元 Gauss 消去法与 Householder 变换法得到的结果相同，二范数都很小。当 n=100 时，列选主元 Gauss 消去法得到的二范数为 0，而 Householder 变换法得到的二范数较大，比前者的误差大得多。当 n=150 时，列选主元 Gauss 消去法得到的二范数还是 0，Householder 变换法得到的二范数也继续增大，误差更大。由此说明，LU 分解列选主元 Gauss 消去法和 QR 分解 Householder 变换法，在不同的条件下误差是不一样的，不仅与 A 矩阵及其阶数 n 有关，也与 b 向量有关。根据实验结果，总体上列选主元 Gauss 消去法要好一些。对线性代数方程组求解，理论方法可以理解，但是在转化为程序代码的时候，花了很大工夫，主要是用 C++的矩阵、向量存储和运算较为复杂，编写程序时几乎都是对矩阵每个元素的变化逐一分析，在这个过程中思路一定要清晰，一旦错误就必须再次分析找到问题所在。最后在实验验证过程中，由于 H、G 矩阵的阶数都很大，不能手动填值，所以要根据规律写相应的矩阵生成函数，简化操作。最后一组实验的理论值没有给出，我通过 MATLAB 编写程序，求得分数形式的结果，作为精确理论值，然后再在 C++程序里求得相应的二范数。

资料库

北邮数值与符号计算实验线性代数方程组求解.docx

相关推荐

课程资源

热门标签

最新资料