《数值计算方法》实验报告
25
X([2,7,12])=1;
d=[1 2 3 4 1 3 2 4]; for i=0:1:100
r1=[cos(i*0.15/100) -sin(i*0.15/100) 0;0 1 0;-sin(i*0.15/100) 0 cos(i*0.15/100)]; U=X*r1';
plot3(U(d,1),U(d,2),U(d,3)); drawnow end
for i=0:1:100
r2=[cos(-i*1.5/100) -sin(-i*0.15/100) 0;sin(-i*0.15/100) cos(-i*0.15/100) 0;0 0 1]; W=U*r2';
plot3(W(d,1),W(d,2),W(d,3)); drawnow end
for i=0:1:100
r3=[1 0 0;0 cos(i*2.7/100) -sin(i*2.7/100) ;0 sin(i*2.7/100) cos(i*2.7/100)]; T=W*r3';
plot3(T(d,1),T(d,2),T(d,3)); drawnow end
subplot(2,2,1);
plot3(X(d,1),X(d,2),X(d,3)); subplot(2,2,2);
plot3(U(d,1),U(d,2),U(d,3)); subplot(2,2,3);
plot3(W(d,1),W(d,2),W(d,3)); subplot(2,2,4);
plot3(T(d,1),T(d,2),T(d,3)); view(3); rotate3d; xlabel('x') ylabel('y') zlabel('z')
P108.1:
function X=gseid2(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
《数值计算方法》实验报告 26
% -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)*P(2))/A(1,1); elseif
X(N)=(B(N)-A(N,N-1)*(X(N-1))')/A(N,N); else
%X contains the kth approximations and P the (k-1)st X(j)=(B(j)-A(j,j-1)*X(j-1)'-A(j,j+1)*P(j+1))/A(j,j); end end
err=abs(norm(X'-P));
relerr=err/(norm(X)+eps); P=X';
if(err %上三角变换 function X=uptrbk(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 and the temporary storage matrix C [N N]=size(A); X=zeros(N,1); C=zeros(1,N+1); %Form the augmented matrix:Aug=[A|B] Aug=[A B]; for p=1:N-1 %Partial proviting for column p [Y,j]=max(abs(Aug(p:N,p))); 《数值计算方法》实验报告 27 %interchange row p and j C=Aug(p,:); Aug(p,:)=Aug(j+p-1,:); Aug(j+p-1,:)=C; if Aug(p,p)==0 'A was singular. No unique solution' break end %Elimination process for cloumn p for k=p+1:N m=Aug(k,p)/Aug(p,p); Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1); end end ock substitustion on [U|Y] using program 3.1 X=backsub(Aug(1:N,1:N),Aug(1:N,N+1)); %回代 function X=backsub(A,B) %input--A is an n*n upper-triangular nonsinggular matrix % -B is an n*1 matrix %output-X is the solution to the linear system AX=B %Find the dimension of B and initialize X n=length(B); X=zeros(n,1); X(n)=B(n)/A(n,n); for k=n-1:-1:1 X(k)=(B(k)-A(k,k+1:n)*X(k+1:n))/A(k,k); end (a)输入： A=[2 0 1;3 2 5;1 -1 0]; B=[1;0;0]; X1=uptrbk(A,B); B=[0;1;0]; X2=uptrbk(A,B); B=[0;0;1]; X3=uptrbk(A,B); [X1 X2 X3] inv(A) (b)输入： A=[16 -120 240 -140;-120 1200 -2700 1680;240 -2700 6480 -4200;-140 1680 -4200 2800]; 《数值计算方法》实验报告 28 B=[1;0;0;0]; X1=uptrbk(A,B); B=[0;1;0;0]; X2=uptrbk(A,B); B=[0;0;1;0]; X3=uptrbk(A,B); B=[0;0;0;1]; X4=uptrbk(A,B); [X1 X2 X3 X4] inv(A) p120.3 function m=newlufact(A,E) %input--A is an N*N nonsingular matrix % -E is an N*N eyes %output-m is an N*N invmatrix of A %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 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;