MATLAB课内实验2程序设计参考.docxVIP

《MATLAB与数学建模》课内实验程序设计

课内实验2

2.1高斯-赛德尔迭代

MATLAB解线性方程组Ax=b，（1）x=A\b；（2）x=inv(A)*b。

一、最简形式1(while循环)：

function[x,k]=GaussSeidelzuijian(A,b,x0,eps,M)

L=tril(A);

U=triu(A,1);

B=-inv(L)*U;%B为迭代矩阵

g=inv(L)*b;%g为右端项

k=0;

tol=1;

whiletol=eps

x=B*x0+g;

k=k+1;%迭代步数。

tol=norm(x-x0);%前后两步迭代结果的误差

x0=x;

if(k=M)

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

return;

end

end

二、最简形式2(for循环)：

function[x,k]=GaussSeidelzuijian2(A,b,x0,eps,M)

L=tril(A);

U=triu(A,1);

B=-inv(L)*U;%B为迭代矩阵

g=inv(L)*b;%g为右端项

k=0;

tol=1;

fork=1:M

x=B*x0+g;

tol=norm(x-x0);%前后两步迭代结果的误差

x0=x;

iftoleps

disp(‘迭代收敛，得到满足精度的解’)

break;

end

end

ifk==M

disp(‘迭代到最大步数M，没有得到满足精度的解’)

end

主程序1（脚本文件fczqj1.m或函数文件）

clear;clc;

A=[-4,1,1,1;1,-4,1,1;1,1,-4,1;1,1,1,-4];b=[1;1;1;1];x0=[0;0;0;0];eps=10^-6;M=1000;

[x,k]=GaussSeidelzuijian(A,b,x0,eps,M)

%[x,k]=GaussSeidelzuijian2(A,b,x0,eps,M)

主程序2（增加注释，完善功能）fczqj2.m

clc;clear;

三、使用稀疏存储，提高计算速度。（求解大型稀疏方程组）

function[x]=xsGs(A,b,x0,e,N)

%功能：用向量稀疏存储形式的G-S迭代法解线性方程组Ax=b

%x0表示迭代初值

%e表示计算精度

%N表示最大迭代次数

%x表示方程组的解向量

%A为系数矩阵

%b为常数向量（方程组右端项）

%

n=length(b);

ifnargin5

N=500;

end

ifnargin4

e=1e-4;

end

ifnargin3

x0=zeros(n,1);

end

x0=sparse(x0);A=sparse(A);b=sparse(b);%采用稀疏存储

x=x0;x0=x+2*e;x0=sparse(x0);

k=0;A1=tril(A);iA1=inv(A1);U=A-A1;

Bgs=-iA1*U;f=iA1*b;

whilenorm(x0-x,inf)ekN,

k=k+1;

x0=x;

x=Bgs*x0+f;

end

x=full(x);

ifk==N,

warning(‘已达迭代次数上限’);

end

注：姜健飞等，数值分析及其MATLAB实验，科学出版社，2004.06，p65。

四、来自于网络：

function[x,k]=GaussSeidel(A,b,x0,eps,M)

%高斯赛德尔迭代法求方程组的解（矩阵公式求解）

%A为方程组的系数矩阵；b为方程组的右端项

%x为线性方程组的解了；x0为迭代初值

%eps为误差限；M为迭代的最大次数

ifnargin==3

eps=1.0e-6;%默认精度

M=10000;%参数不足时默认后两个条件

elseifnargin==4

M=10000;%参数的默认值

elseifnargin3

error(参数不足);

return

end

[n,m]=size(A);

nb=length(b);

%当方程组行与列的维数不相等时，停止计算，并输出出错信息

ifn~=m

error(矩阵A行数和列数必须相等!);

return;

end

%当方程组与右端项的维数不匹配时，停止计算，并输出出错信息

ifn~=nb

error(矩阵A的行数必须和b的长度相等!);

return;

end

L=zeros(n,n)