作为研究生的入门课，数值计算的大作业算是所有研究生开学的重要编程作业。

我把矩阵直接求解的LU分解法与Chloesky分别在MATLAB和C语言中编程实现。具体的程序详细标注后放在文章最后了，需要的同学自取。

LU分解数学原理：

将系数矩阵A转变成等价两个矩阵L和U的乘积，其中L和U分别是单位下三角矩阵和上三角矩阵。当A的所有顺序主子式都不为0时，矩阵A可以分解为A=LU。其中L是下三角矩阵，U是上三角矩阵。这时方程组Ax=b可以写成LUX=b.

令UX=Y,则有LY=b。由此，原方程组求解变为，求解系数矩阵为三角阵的方程，很容易实现。具体计算步骤如下：

计算A的LU分解
向前回代求解下三角方程LY=b
向后回代求解上三角方程UX=Y

下文为作业详解

题目：考虑一个一端固定，另一端自由的水平悬杆，其受力的离散模型产生了一个五对角线性代数方程组Ax=b，其中 n维向量b为杆上的已知荷载，假定其为均匀，即b的每个分量为1，n维向量x表示杆的偏移。完成以下工作

1.使用LU分解，分别求解n=100和n=1000时方程组的数值解，并计算相应最大残差，说明求得的解是否准确？并从系数矩阵条件数出发给予简单解释。

具体C语言与Matlab程序在附录中，此处只展示结果与结论

当n=100时，部分解的运算结果如下

列 1 至 14

0.0001 0.0007 0.0019 0.0034 0.0055 0.0080 0.0109 0.0143 0.0181 0.0224 0.0270 0.0321 0.0375 0.0433

列 15 至 28

0.0495 0.0561 0.0630 0.0702 0.0778 0.0858 0.0940 0.1026 0.1115 0.1207 0.1302 0.1399 0.1500 0.1603

列 29 至 42

0.1708 0.1817 0.1927 0.2041 0.2156 0.2274 0.2394 0.2516 0.2640 0.2766 0.2895 0.3025 0.3157 0.3290

.............

列 85 至 98

1.0025 1.0193 1.0360 1.0528 1.0695 1.0863 1.1031 1.1199 1.1366 1.1534 1.1702 1.1870 1.2038 1.2206

列 99 至 100

1.2374 1.2542

相应的最大残差residual=7.4506*10^(-9)

系数矩阵A的条件数为2.0067*10^8(远大于1）

当n=1000时，部分解的运算结果如下

列 1 至 14

0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001 0.0002 0.0002 0.0003 0.0003 0.0004 0.0005

列 15 至 28

0.0005 0.0006 0.0007 0.0008 0.0009 0.0010 0.0011 0.0012 0.0013 0.0014 0.0015 0.0016 0.0018 0.0019

...........

列 967 至 980

1.1954 1.1970 1.1987 1.2004 1.2020 1.2037 1.2054 1.2071 1.2087 1.2104 1.2121 1.2137 1.2154 1.2171

列 981 至 994

1.2187 1.2204 1.2221 1.2237 1.2254 1.2271 1.2287 1.2304 1.2321 1.2337 1.2354 1.2371 1.2387 1.2404

列 995 至 1,000

1.2421 1.2437 1.2454 1.2471 1.2488 1.2504

相应的最大残差residual=1.0681*10^(-4)

系数矩阵A的条件数为2.0007*10^12(远大于1)

结论：根据上述两个不同的系数矩阵LU分解的解，当系数矩阵从100维放大到1000维时，最大残差明显放大了10^5。这是由于数值分析中，一个问题的条件数是该数量在数值计算中的容易程度的衡量，也就是该问题的适定性。一个低条件数的问题称为良态的，而高条件数的问题称为病态（或者说非良态）的。条件数就是一个矩阵的一种范数乘以此逆矩阵的范数。（本次实验计算条件数为系数矩阵的无穷范数）因此条件数越大，矩阵A或常数项b引起方程组的解变化越巨大。所以上述n=100的系数矩阵条件数远小于n=1000时的条件数，故n=100时的数值解更接近实际值，最大残差也更小。

2.由于系数矩阵对称正定，尝试给出矩阵当 n=5和n=10时的Cholesky分解，并求解方程组，计算相应最大残差

当n=5时候，利用Cholesky分解求得的解以及最大残差如下图所示

当n=10时候，利用Cholesky分解求得的解以及最大残差如下图所示

附录1：Matlab程序直接法求解线性方程组

% 线性方程组数值求解
% 实验一直接法解方程
% 产生n阶带状稀疏矩阵A，b
n=10;
b=ones(1,n)';
%生成所有对角元素的向量(一共有5个长度为n的对角向量，分别命名为1，2，3，4，5，其中主对角线向量为diagA_3)
diagA_1=zeros(1,n)';
diagA_2=zeros(1,n)';
diagA_3=zeros(1,n)';
diagA_4=zeros(1,n)';
diagA_5=zeros(1,n)';
% 对角线-2赋值for i=1:n-2diagA_1(i)=1;end
% 对角线-1赋值
diagA_2(n-1)=-2;
for i=1:n-2diagA_2(i)=-4;
end
% 主对角线赋值
diagA_3(1)=9;
diagA_3(n-1)=5;
diagA_3(n)=1;
for i=1:n-3diagA_3(i+1)=6;
end
% 对角线+1赋值
diagA_4(n)=-2;
for i=1:n-2diagA_4(i+1)=-4;
end% 对角线+2赋值
for i=1:n-2diagA_5(i+2)=1;
end
% 将5个对角向量生成一个过渡矩阵B，以及带状位置构成向量d
B=[diagA_1,diagA_2,diagA_3,diagA_4,diagA_5];
d=[-2,-1,0,1,2]';
A=spdiags(B,d,n,n);
A=full(A);% LU分解求方程解[l,u]=lu(A);x=u\(l\b);% 计算A矩阵在各种范数条件下的条件数% A的1-范数条件数目cond_num1=cond(A,1);% A的2-范数条件数目cond_num2=cond(A,2);% A的-范数条件数目cond_num3=cond(A,inf);% % Cholesky分解法求解方程
% ch=chol(A);
% x=ch\(ch'\b);
% 计算相应的最大残差C=b-A*x;residual=norm(C,inf);
%  展示最后计算的数值解x,对应的最大残差和系数矩阵的条件数
disp(x');
disp(residual);
disp(cond_num3);

附录2：C语言程序LU按行列分别求解方程组（只尝试LU分解部分）

逐步进行LU分解
for(i=0; i<n-1; i++)
{/*对“L列”进行计算*/for(j=i+1; j<n; j++)
{for(k=0,sum=0; k<n; k++){if(k != i) sum += l[j][k]*u[k][i];}l[j][i] = (float)((a[j][i]-sum)/u[i][i]);}/*对“U行”进行计算*/
for(j=i+1; j<n; j++){for(k=0,sum=0; k<n; k++){if(k != i+1) sum += l[i+1][k]*u[k][j];}u[i+1][j] = (float)((a[i+1][j]-sum));}}
/*输出矩阵l*/
printf("矩阵L：\n");
for(i=0; i<n; i++){for(j=0; j<n; j++){printf("%0.3f  ", l[i][j]);}printf("\n");}
/*输出矩阵u*/
printf("\n矩阵U：\n");
for(i=0; i<n; i++)
{for(j=0; j<n; j++){printf("%0.3f  ", u[i][j]);}printf("\n");
}

数值计算大作业：线性方程组求解（LU分解与Chloesky分解程序在Matlab与C的实现）相关推荐

数值计算大作业：Jacobi与Gauss -Seidel迭代求解线性方程组（Matlab实现）
作为研究生的入门课,数值计算的大作业算是所有研究生开学的重要编程作业. 我把Jacobi与Gauss -Seidel迭代求解线性方程组的数值计算作业在MATLAB中编程实现.具体的程序详细标注后放在文 ...
李保滨矩阵分析大作业2022：LU、QR、URV分解、Householder、Givens变换的程序实现
介绍本文为2022年秋季学期国科大李保滨老师的矩阵分析与应用课程大作业实现,编程语言使用python 具体作业要求: 完成课堂上讲的关于矩阵分解的LU.QR(Gram-Schmidt).正交规约(H ...
数值计算大作业：最小二乘法拟合（Matlab实现）
作为研究生的入门课,数值计算的大作业算是所有研究生开学的重要编程作业. 我把最小二乘算法在MATLAB中整合成了一个M函数文件least square fitting.m,直线拟合函数lsf_line ...
数值计算大作业：非线性方程求根（二分法、牛顿法、弦截法在Matlab实现）
作为研究生的入门课,数值计算的大作业算是所有研究生开学的重要编程作业. 我把二分法.牛顿法.弦截法求解非线性方程求根的数值计算作业在MATLAB中编程实现.具体的程序详细标注后放在文章附录了,算法数学 ...
数值计算大作业：数值积分（梯形、辛普森与龙贝格方法在Matlab实现）
作为研究生的入门课,数值计算的大作业算是所有研究生开学的重要编程作业. 针对数值积分的编程,我把梯形.辛普森与龙贝格方法在MATLAB中编写的程序放在文章最后了,需要的同学自取. PS:附录中的程序并 ...
java大作业国际比赛奖牌榜,Java中集合的程序练习
一.运用所学集合的知识,写出一个奥运会奖牌榜排序程序,具体要求如下: 1.每个国家都分别拥有金牌银牌铜牌属性. 2.对各个国家实现奖牌榜排序,排序规则为:先比较金牌,如果金牌数相同则比较银牌,如果银牌 ...
2022春哈工大计算机系统大作业——hello的程序人生
计算机系统大作业题目程序人生-Hello's P2P 专业计算学部学号班级学生指导教师计算机科学与技术学院 2021年5月摘 ...
计算机系统大作业————程序人生
计算机系统大作业题目程序人生-Hello's P2P 专业计算机学号 *********** 班级 ******* 学生 *** 指导教师郑贵滨计算机科学与技术学院 20 ...
python面向对象课程大作业定义一个描述学生基本情况的类，数据成员至少包括 “姓名、性别、学号、年级、所在院系、面向对象的考试日期”
python面向对象课程大作业按下列要求编写一个完整的程序: 定义一个描述学生基本情况的类,数据成员至少包括"姓名.性别.学号.年级.所在院系.面向对象的考试日期",成员函数至少 ...

数值计算大作业：线性方程组求解（LU分解与Chloesky分解程序在Matlab与C的实现）