数值分析上机实验——解线性方程组(2)

end

②控制台输入代码： >> a=[2 2 2 2 2]; >> b=[-1 -1 -1 -1]; >> c=[-1 -1 -1 -1]; >> f=[1;0;0;0;0]; >> x=ZG_SDJ(a,b,c,f) 3.

①改进的平方根法（文件GJPFG.m） function GJPFG(A,b) n=length(b);% n?aáD??￡? % LDL'·??a￡? d(1)=A(1,1); for i=2:n

for j=1:i-1 sum1=0; for k=1:j-1

sum1=sum1+t(i,k)*l(j,k); end

t(i,j)=A(i,j)-sum1; l(i,j)=t(i,j)/d(j); end sum2=0; for k=1:i-1

sum2=sum2+t(i,k)*l(i,k); end

d(i)=A(i,i)-sum2; end for i=1:n l(i,i)=1; end

disp('μ￥????èy?????óL?a￡o'); %?a3?μ￥????èy?????óL￡? l

disp('???????óD?a￡o'); %?a3????????óD￡? d

%óéLDL'x=b?ó?ax￡? %óéLy=b￡??óy￡?

%óéL'x=inv￡¨D￡?y￡??ó?ax￡? y(1)=b(1); for i=2:n sum3=0; for k=1:i-1

sum3=sum3+l(i,k)*y(k);

end

y(i)=b(i)-sum3; end

x(n)=y(n)/d(n); for i=n-1:-1:1 sum4=0; for k=i+1:n

sum4=sum4+l(k,i)*x(k); end

x(i)=(y(i)/d(i))-sum4; end

disp('óéLy=b?ó?ayμ?￡o'); y

disp('Ax=bμ??ax?a￡o'); x

②控制台输入代码：

>> A=[2 -1 1;-1 -2 3;1 3 1]; >> b=[4;5;6]; >> GJPFG(A,b) 4.

①SOR方法（文件SOR_1.m） function SOR_1(A,b,x0,x_a,w) %x_a?a??è·?a

if(w<=0 || w>=2)

error('2?êy·??§′í?ó'); return; end

eps=5.0e-6;

D=diag(diag(A)); %?óAμ????????ó L=-tril(A,-1); %?óAμ???èy???ó U=-triu(A,1); %?óAμ?é?èy???ó B=inv(D-L*w)*((1-w)*D+w*U); f=w*inv((D-L*w))*b; x=B*x0+f;

n=1; %μü′ú′?êy

while norm(x_a-x)>=eps x0=x; x =B*x0+f; n=n+1; if(n>=200)

disp('Warning: μü′ú′?êyì??à￡??é?ü2?ê?á2￡?');

return; end end

disp('Ax=bμ??a?a￡o'); x

disp('μü′ú′?êy?a￡o'); n

②控制台输入代码：

>> A=[4 -1 0;-1 4 -1;0 -1 4]; >> b=[1;4;-3]; >> x0=[0;0;0];

>> x_a=[0.5;1;-0.5]; >> w=1.03;

>> SOR_1(A,b,x0,x_a,w) >> w=1;

>> SOR_1(A,b,x0,x_a,w) >> w=1.1;

>> SOR_1(A,b,x0,x_a,w) 5．

①SOR方法（文件SOR_2.m） function SOR_2(A,b,x0,w,eps) if(w<=0 || w>=2)

error('2?êy·??§′í?ó'); return; end

D=diag(diag(A)); %?óAμ????????ó L=-tril(A,-1); %?óAμ???èy???ó U=-triu(A,1); %?óAμ?é?èy???ó B=inv(D-L*w)*((1-w)*D+w*U); f=w*inv((D-L*w))*b; x=B*x0+f;

n=1; %μü′ú′?êy

while norm(x-x0)>=eps x0=x; x =B*x0+f; n=n+1; if(n>=200)

disp('Warning: μü′ú′?êyì??à￡??é?ü2?ê?á2￡?'); return; end end

disp('Ax=bμ??a?a￡o'); x

disp('μü′ú′?êy?a￡o'); n

②控制台输入代码：

>> A=[5 2 1;-1 4 2;2 -3 10]; >> b=[-12;20;3]; >> x0=[0;0;0]; >> w=0.9;

>> eps=10e-4;

>> SOR_2(A,b,x0,w,eps)

6．此题为证明题，无程序代码。 7.

① 雅克比迭代法（文件Jocabi.m） function x=Jocabi(n) A=hilb(n);

x_a=ones(n,1); b=A*x_a; eps=1e-4;

n=length(b); N=50;

x=zeros(n,1); y=zeros(n,1);

for k=1:N sum=0; for i=1:n

y(i)=(b(i)-A(i,1:n)*x(1:n)+A(i,i)*x(i))/A(i,i); end

for i=1:n

sum=sum+(y(i)-x(i))^2; end

if sqrt(sum)

for i=1:n x(i)=y(i); end end

end

② SOR方法（文件SOR_3.m） function SOR_3(n,w) %x_a?a??è·?a

if(w<=0 || w>=2)

error('2?êy·??§′í?ó'); return; end

x0=zeros(n,1); A=hilb(n);

x_a=ones(n,1); b=A*x_a; eps=1e-4;

D=diag(diag(A)); %?óAμ????????ó L=-tril(A,-1); %?óAμ???èy???ó U=-triu(A,1); %?óAμ?é?èy???ó B=inv(D-L*w)*((1-w)*D+w*U); f=w*inv((D-L*w))*b; x=B*x0+f;

n=1; %μü′ú′?êy

while norm(x-x0)>=eps x0=x; x =B*x0+f; n=n+1; if(n>=2000)

disp('Warning: μü′ú′?êyì??à￡??é?ü2?ê?á2￡?'); return; end end

disp('Hx=bμ??a?a￡o'); x

disp('μü′ú′?êy?a￡o'); n

③控制台输入代码： >> x=Jocabi(6) >> x=Jocabi(8) >> x=Jocabi(10) >> SOR_3(6,1) >> SOR_3(6,1.25) >> SOR_3(6,1.5)

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库数值分析上机实验——解线性方程组(2)在线全文阅读。