资源描述:
《计算方法实验报告习题2(浙大版)》由会员上传分享,免费在线阅读,更多相关内容在教育资源-天天文库。
1、专业:电气工程及其自动化姓名:李X计算方法实验报告实验名称:实验2列主元素消去法解方程组1引言工程实际问题中,线型方程的系数矩阵一般为低阶稠密矩阵和大型稀疏矩阵。用高斯消去法解Ax=b时,可能出现很小,用作除数会导致中间结果矩阵元素数量级严重增长和舍入误差的扩散,使结果不可靠;采用选主元素的三角分解法可以避免此类问题。高斯消去法的消去过程,实质上是将A分解为两个三角矩阵的乘积A=LU,并求解Ly=b的过程。回带过程就是求解上三角方程组Ux=y。所以在实际的运算中,矩阵L和U可以直接计算出,而不需要任何中间步骤,从而在计算过程中
2、将高斯消去法的步骤进行了进一步的简略,大大提高了运算速度,这就是三角分解法。采用选主元的方式与列主元高斯消去法一样,也是为了避免除数过小,从而保证了计算的精确度。2实验目的和要求通过列主元素消去法求解线性方程组,实现PA=LU。要求计算解x,L,U,整形数组IP(i),(i=1,2,…,)(记录主行信息)。3算法原理与流程图(1)原理将A分解为两个三角矩阵的乘积A=LU。对方程组的增广矩阵经过k-1步分解后,可变成如下形式:第k步分解,为了避免用绝对值很小的数作除数,引进量,于是有=。如果,则将矩阵的第t行与第k行元素互换,将
3、(i,j)位置的新元素仍记为或,然后再做第k步分解,这时(2)流程图1程序代码及注释%运用列主元素消去法求解方程组并实现PA=LUfunction[x,L,U,P,IP]=lzylu(A,b)%定义函数[m,n]=size(A);%计算系数矩阵A的行列数ifm~=nerror('系数矩阵不是方阵');return;endifdet(A)==0%计算矩阵A的行列式,若为零则A是奇异矩阵error('矩阵不能被三角形分解')endu=A;v=b;q=eye(m);s=q;%定义单位矩阵forj=1:m-1%按照列循环t=abs(u
4、(j,j));mj=j;fori=j:mifabs(u(i,j))>tmj=i;endend%比较大小,选列主元素ifu(mj,j)==0error('系数矩阵奇异');returnendB=u(j,:);c=v(j);u(j,:)=u(mj,,:);v(j)=v(mj);u(mj,:)=B;v(mj)=c;%交换两行,将列主元素所在行移到第j行IP(j)=mj;%记录互换,表示第j行与第IP(j)行互换ts=s(j,:);s(j,:)=s(mj,:);s(mj,:)=ts;%置换矩阵进行相应行交换u(j+1:m,j)=u(j
5、+1:m,j)/u(j,j);%求出L矩阵第j列第j+1到第m行元素u(j+1:m,j+1:m)=u(j+1:m,j+1:m)-u(j+1:m,j)*u(j,j+1:m);%解出u矩阵做完第j次消元后的矩阵v(j+1:m)=v(j+1:m)-v(j)*u(j+1:m,j);%v也要对应系数进行运算,相比杜利特尔(Doolittle)分解法减少了l*y=b,y=u*x先要求解y再求x的过程endP=s;L=tril(u,-1)+q;%对u取下三角,其余部分全部为0,得到L矩阵U=triu(u);%对u取上三角,其余部分全部为0,
6、得到U矩阵x(m)=v(m)/u(m,m);fori=m-1:-1:1sum=0;fork=i+1:msum=sum+U(i,k)*x(k);endx(i)=(v(i)-sum)/u(i,i);end%利用U*x=v,求解x得到行向量x=x';%转置为列向量1算例分析解方程组>>A=[1111111;2111111;3211111;4321111;5432111;6543211;7654321];>>b=[781013172228]';>>[x,L,U,P,IP]=lzylu(A,b)x=1.00001.00001.00001
7、.00001.00001.00001.0000L=1.00000000000.28571.0000000000.42860.80001.000000000.57140.60000.75001.00000000.71430.40000.50000.66671.0000000.1429-0.2000-0.2500-0.3333-0.50001.000000.85710.20000.25000.33330.5000-1.00000.0000U=7.00006.00005.00004.00003.00002.00001.00000-0
8、.7143-0.4286-0.14290.14290.42860.714300-0.8000-0.6000-0.4000-0.20000.0000000-0.7500-0.5000-0.25000.00000000-0.6667-0.3333-0.0000000000.500