2013级工科硕士研究生

《矩阵与数值分析》课程数值实验题目

一、设SN??j?2N106，分别编制从小到大和从大到小的顺序程序分别计算 j2?1S10000,S1000000

并指出两种方法计算结果的有效位数。 Matlab程序如下：

function [si,sd]=S(N) format long; si=0;sd=0; for j=N:-1:2

si=1.0e6/(j^2-1)+si; end

for j=2:N

sd=1.0e6/(j^2-1)+sd; end end

在matlab命令窗口中输入：[si,sd]=S(10000) 运行结果：si =7.499000049995000e+005

sd =7.499000049994994e+005

在matlab命令窗口中输入：[si,sd]=S(1000000) 运行结果：si =7.499990000005000e+005

sd =7.499990000005200e+0051

结果分析：si为从大到小的顺序求和的值，sd为从小到大的顺序求和的值。当N分别为10000和1000000时，si分别为7.499000049995000e+005和7.499990000005000e+005，可以看出这两个数的有效值均为13位；而sd分别为7.499000049994994e+005和7.499990000005200e+005，这两个数的有效值均为16位。这就出现了我们在矩阵理论课上所学的“大数吃小数”的问题。为了使结果更为精确我们必须避免在四则运算中出现“大数吃小数”的情况，应该按从小到大的顺序进行求和。

二、解线性方程组

1．分别利用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组Ax?b，其中常向量为?n?1?维随机生成的列向量，系数矩阵A具有如下形式

2?Tn?1?2In?1?In?1????In?1?， A????In?1????In?1Tn?1?2In?1???2?1????1其中T??迭代法计算停止的条件?为n?1阶矩阵，In?1为n?1阶单位矩阵，n?1??1????12??为：xk?1?xk?10，给出n?10,100,1000时的不同迭代步数． ?102求解系数矩阵A和向量b的Matlab程序如下：

function [A,b,x0]=jz(n)

for i=1:n-1 %此处n赋值即可，如n=100 for j=1:n-1 if(i==j)

T(i,j)=2; end

if(i==(j+1)) T(i,j)=-1; end

if(i==(j-1)) T(i,j)=-1; end end end

I=eye(n-1); k=1;

for mm=1:(n-1)

A(k:(k+n-2),k:(k+n-2))=T+2*I; k=k+n-1; end k=1;

for xx=1:(n-2)

A(k:(k+n-2),(k+n-1):(k+2*n-3))=-I; A((k+n-1):(k+2*n-3),k:(k+n-2))=-I; k=k+n-1; end

b=randn((n-1)^2,1); x0=zeros((n-1)^2,1);

Jacobi迭代法Matlab程序如下：

function n=jacobi(A,b,x0) eps= 1.0e-10;

M = 10000;

D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 B=D\\(L+U); f=D\\b; x=B*x0+f;

n=1; %迭代次数 while norm(x-x0)>=eps x0=x;

x=B*x0+f; n=n+1; if(n>=M)

disp('Warning: 迭代次数太多，可能不收敛！'); return; end end

Gauss-Seidel迭代法Matlab程序如下：

function n=gauseidel(A,b,x0) eps= 1.0e-10; M = 10000;

D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 G=(D-L)\\U; f=(D-L)\\b; x=G*x0+f;

n=1; %迭代次数 while norm(x-x0)>=eps x0=x;

x=G*x0+f; n=n+1; if(n>=M)

disp('Warning: 迭代次数太多，可能不收敛！'); return; end end

在matlab命令窗口中输入：

[A,b,x0]=jz(10); J10=jacobi(A,b,x0) G10=gauseidel(A,b,x0) [A,b,x0]=jz(20); J20=jacobi(A,b,x0) G20=gauseidel(A,b,x0)