南京邮电大学 数值代数实验

南京邮电大学 数值代数实验

ID:1243607

大小:1.44 MB

页数:21页

时间:2017-11-09

南京邮电大学 数值代数实验_第1页
南京邮电大学 数值代数实验_第2页
南京邮电大学 数值代数实验_第3页
南京邮电大学 数值代数实验_第4页
南京邮电大学 数值代数实验_第5页
资源描述:

《南京邮电大学 数值代数实验》由会员上传分享,免费在线阅读,更多相关内容在行业资料-天天文库

1、数值代数实验数值线性代数实验一一、实验名称:矩阵的LU分解.二、实验目的:用不选主元的LU分解和列主元LU分解求解线性方程组Ax=b,并比较这两种方法.三、实验内容与要求(1)用所熟悉的计算机语言将不选主元和列主元LU分解编成通用的子程序,然后用编写的程序求解下面的84阶方程组将计算结果与方程组的精确解进行比较,并就此谈谈你对Gauss消去法的看法.(2)写出追赶法求解三对角方程组的过程,并编写程序求该实验中的方程组Gauss消去法:用消去法解方程组的基本思想是用逐次消去未知数的方法把原来方程组Ax=b化为与其等价的三角方程组,而求

2、解三角方程组就容易了。换句话说,上述过程就是用行的初等变换将原方程组系数矩阵化为简单形式,从而将求解原方程组的问题转化为求解简单方程组的问题。利用Gauss消去法对线性方程组Ax=b进行求解。用MATLAB建立m文件DelGauss.m,程序如下:functionx=DelGauss(a,b)[n,m]=size(a);nb=length(b);det=1; x=zeros(n,1);fork=1:n-1fori=k+1:nifa(k,k)==0returnendm=a(i,k)/a(k,k);forj=k+1:na(i,j)=a(

3、i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k); enddet=det*a(n,n);fork=n:-1:1forj=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);End在matlab中输入如下:结果如下:方程组的精确解为x1=x2=…=x84=1.0000,与Gauss消去法求得的解差距很大,所得结果不够准确,计算简单但其消元过程有时不能进行到底而使求解出现解失真的情况。数值线性代数实验二一、实验名称:实对称正定矩阵的的Chol

4、esky分解.二、实验目的:用平方根法和改进的平方根方法求解线性方程组Ax=b.三、实验内容与要求用所熟悉的计算机语言将Cholesky分解和改进的Cholesky分解编成通用的子程序,然后用编写的程序求解对称正定方程组Ax=b,其中(1)b随机的选取,系数矩阵为100阶矩阵(2)系数矩阵为40阶Hilbert矩阵,即系数矩阵A的第i行第j列元素为,向量b的第i个分量为(3)用实验一的程序求解这两个方程组,并比较所有的计算结果,然后评价各个方法的优劣。平方根法:平方根法就是利用对称正定矩阵的三角分解而得到的求解对称正定方程组的一种有

5、效方法。平方根法递推公式可以证明对于对称正定矩阵A,可以唯一地分解成A=LLT,其中L是非奇异下三角形矩阵。模型二:利用平方根法对线性方程组Ax=b进行求解。用MATLAB建立m文件pingfg.m,程序如下:function[x]=pingfg(A,b)%Cholesky分解[n,n]=size(A);L=zeros(n,n);%实际上不用为L申请空间,使用A即可L(1,1)=sqrt(A(1,1));fork=2:nL(k,1)=A(k,1)/L(1,1);endfork=2:n-1L(k,k)=sqrt(A(k,k)-sum(

6、L(k,1:k-1).^2));fori=k+1:nL(i,k)=(A(i,k)-sum(L(i,1:k-1).*L(k,1:k-1)))/L(k,k);endendL(n,n)=sqrt(A(n,n)-sum(L(n,1:n-1).^2));%解下三角方程组Ly=by=zeros(n,1);fork=1:nj=1:k-1;y(k)=(b(k)-L(k,j)*y(j))/L(k,k);end%解上三角方程组L'x=yx=zeros(n,1);U=L';fork=n:-1:1j=k+1:n;x(k)=(y(k)-U(k,j)*x(j)

7、)/U(k,k);End模型三:利用改进的平方根法对线性方程组Ax=b进行求解。用MATLAB建立m文件ave.m,程序如下:function[x]=ave(A,b,n)%用改进平方根法求解Ax=bL=zeros(n,n);%L为n*n矩阵D=diag(n,0);%D为n*n的主对角矩阵S=L*D;fori=1:n%L的主对角元素均为1L(i,i)=1;endfori=1:nforj=1:n%验证A是否为对称正定矩阵if(eig(A)<=0)

8、(A(i,j)~=A(j,i))%A的特征值小于0或A非对称时,输出wrongdisp('

9、wrong');break;endendendHilbert矩阵用MATLAB建立m文件Hil.m,程序如下:functionb=Hil()fork=1:40form=1:40s=0;t=s+1/(k+m-1);s=t;endb(k,

当前文档最多预览五页,下载文档查看全文

此文档下载收益归作者所有

当前文档最多预览五页,下载文档查看全文
温馨提示:
1. 部分包含数学公式或PPT动画的文件,查看预览时可能会显示错乱或异常,文件下载后无此问题,请放心下载。
2. 本文档由用户上传,版权归属用户,天天文库负责整理代发布。如果您对本文档版权有争议请及时联系客服。
3. 下载前请仔细阅读文档内容,确认文档内容符合您的需求后进行下载,若出现内容与标题不符可向本站投诉处理。
4. 下载文档时可能由于网络波动等原因无法下载或下载错误,付费完成后未能成功下载的用户请联系客服处理。