《数值计算方法》实验报告

29

A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N); end end

%solve for Y for J=1:N

Y(1)=E(t(1),J); for k=2:N

Y(k)=E(t(k),J)-A(k,1:k-1)*Y(1:k-1); end

%solve for X

X(N)=Y(N)/A(N,N); for k=N-1:-1:1

X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k); end

m(1:N,J)=X(1:N); end

P120.4

function X=lufact(A,B)

%input--A is an N*N nonsingular matrix % -B is an N*1 matrix

%output-X is an N*1 matrix containing the solution to AX=B

%initialize X, Y, the temporary storage matrix C, and the row %permutation information matrix t [N,N]=size(A); X=zeros(N,1); Y=zeros(N,1); C=zeros(1,N); t=1:N;

for p=1:N-1

%find the privot row for column p [max1,j]=max(abs(A(p:N,p))); %interchange row p and j C=A(p,:);

A(p,:)=A(j+p-1,:); A(j+p-1,:)=C; d=t(p);

t(p)=t(j+p-1); t(j+p-1)=d; if A(p,p)==0

'A is singular. No unique solution' break

《数值计算方法》实验报告 30

end

êlculate multiplier and place in subdiagonal protion of A for k=p+1:N

mult=A(k,p)/A(p,p); A(k,p)=mult;

A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N); end end

%solve for Y Y(1)=B(t(1)); for k=2:N

Y(k)=B(t(k))-A(k,1:k-1)*Y(1:k-1); end

%solve for X

X(N)=Y(N)/A(N,N); for k=N-1:-1:1

X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k); end (a)

输入：R1=1;R2=1;R3=2;R4=1;R5=2;R6=4;E1=23;E2=29; A=[R1+R3+R4 R3 R4;R3 R2+R3+R5 -R5;R4 -R5 R4+R5+R6] B=[E1;E2;0] X=lufact(A,B) （b） 输入：

R1=1;R2=0.75;R3=1;R4=2;R5=1;R6=4;E1=12;E2=21.5; A=[R1+R3+R4 R3 R4;R3 R2+R3+R5 -R5;R4 -R5 R4+R5+R6] B=[E1;E2;0] X=lufact(A,B) （c）输入：

R1=1;R2=2;R3=4;R4=3;R5=1;R6=5;E1=41;E2=38;

A=[R1+R3+R4 R3 R4;R3 R2+R3+R5 -R5;R4 -R5 R4+R5+R6] B=[E1;E2;0] X=lufact(A,B)

《数值计算方法》实验报告

31

P130.4

function X=gseid1(A,B,P,delta,max1) %input--A is an N*N nonsingular matrix % -B is an N*1 matrix % -P is an N*1 matrix

% -delta is the tolerance for P

% -max1 is the maximum number of iterations

%output-X is an N*1 matrix:the gauss-seidel approximation to the solution

%of AX=B

N=length(B);

for k=1:max1 for j=1:N if j==1

X(1)=(B(1)-A(1,2:3)*P(2:3))/A(1,1); elseif j==2

X(2)=(B(2)-A(2,1)*X(1)'-A(2,3:4)*P(3:4))/A(2,2); elseif j==N-1

X(N-1)=(B(N-1)-A(N-1,N-3:N-2)*X(N-3:N-2)'-A(N-1,N)*P(N))/A(N-1,N-1);

elseif j==N

X(N)=(B(N)-A(N,N-2:N-1)*(X(N-2:N-1))')/A(N,N); else

%X contains the kth approximations and P the (k-1)st

X(j)=(B(j)-A(j,j-2:j-1)*X(j-2:j-1)'-A(j,j+1:j+2)*P(j+1:j+2))/A(j,j);

end end

err=abs(norm(X'-P));

relerr=err/(norm(X)+eps); P=X';

if(err

A=zeros(50);

《数值计算方法》实验报告 32

A(1,1:3)=[12 -2 1]; A(2,1:4)=[-2 12 -2 1]; for j=3:50

A(j,j-2:j+2)=[1 -2 12 -2 1]; end

A(49,47:50)=[1 -2 12 -2]; A(50,48:50)=[1 -2 12]; B=5.*ones(50,1); p=zeros(50,1); d=1e-10;

x=gseid1(A,B,p,d,1000)