数值分析MATLAB算法程序附录

数值分析MATLAB算法程序附录

Chapter 1

%程序1.1 顺序高斯消元法  xmyguass.m
function [x]=xmygauss(A,b,flag) %输入系数矩阵A，右端项b,是否显示中间过程的参数flag，取值为0,1；输出解x
n=length(b);
for k=1:(n-1)
m=A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n);
b(k+1:n)=b(k+1:n)-m*b(k);
A(k+1:n,k)=zeros(n-k,1);
if flag~=0
C=horzcat(A,b)  %horzcat函数
end
end
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

function [x]=xmylguass(A,b,flag)
if nargin<3,flag=0;end
n=length(b);
for k=1:(n-1)
[ufo,ind]=max(abs(A(k:n,k)));
ind=ind+k-1;
if ind>k
A([k,ind],:)=A([ind,k],:);
b([k,ind],:)=b([ind,k],:);
end
a=A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-a*A(k,k+1:n);
b(k+1:n)=b(k+1:n)-a*b(k);
A(k+1:n,k)=zeros(n-k,1);
if flag~=0
C=horzcat(A,b);
disp(C)
end
if A(n,n)==0,
break;
end
end
x=zeros(n,1);
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
x(i)=(x(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);
end

function Z=xmyinv(A)
n=length(A) % n=size(A,2)
I=eye(n,n);
C=horzcat(A,I)
for k=1:n
[ufo,ind]=max(abs(C(k:n,k)));
ind=ind+k-1;
if ind>k
C([k,ind],:)=C([ind,k],:);
end
for i=1:n
if i~=k
C(i,k)=C(i,k)/C(k,k);
C(k,k)=1/C(k,k);
C(i,k+1:2*n)=C(i,k+1:2*n)-C(i,k)*C(k,k+1:2*n);
C(k,k+1:2*n)=C(k,k+1:2*n)*C(k,k);
end
end
disp(C)
end
Z=C(1:n,n+1:2*n)

function [x,L,U]=xmycrout(A,b)
n=length(b);
L=zeros(n,n);
U=eye(n,n);
L(:,1)=A(:,1);
U(1,2:n)=A(1,2:n)/L(1,1);
for k=2:n
L(k:n,k)=A(k:n,k)-L(k:n,1:k-1)*U(1:k-1,k);
U(k,k+1:n)=(A(k,k+1:n)-L(k,1:k-1)*U(1:k-1,k+1:n))/L(k,k);
end
L,U
y=zeros(n,1);
y(1)=b(1)/L(1,1);
for k=2:n
y(k)=(b(k)-L(k,1:k-1)*y(1:k-1))/L(k,k);
end
x=zeros(n,1);
x(n)=y(n);
for k=n-1:-1:1
x(k)=y(k)-U(k,k+1:n)*x(k+1:n);
end

function [x,L,U,P]=xmydoolittle(A,b)
n=length(b);
U=zeros(n,n);
L=eye(n,n);
U(1,:)=A(1,:);L(2:n,1)=A(2:n,1)/U(1,1);
for k=2:n
U(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);
L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k))/U(k,k);
end
y=zeros(n,1);
y(1)=b(1);
for k=2:n
y(k)=b(k)-L(k,1:k-1)*y(1:k-1);
end
x=zeros(n,1);
x(n)=y(n)/U(n,n);
for k=n-1:-1:1
x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);
end

function [x,L,U,P]=xmylu(A,b)
n=length(b);
P=eye(n);
for k=1:n
A(k:n,k)=A(k:n,k)-A(k:n,1:k-1)*A(1:k-1,k);
[ufo,mmax]=max(abs(A(k:n,k)));
mmax=mmax+k-1;
if mmax~=k
A([k mmax],:)=A([mmax k],:);
P([k mmax],:)=P([mmax k],:);
b([k mmax],:)=b([mmax k],:);
end
A(k+1:n,k)=A(k+1:n,k)/A(k,k);
A(k,k+1:n)=A(k,k+1:n)-A(k,1:k-1)*A(1:k-1,k+1:n);
end
L=tril(A,-1)+eye(n,n); % L(k+1:n,k)=A(k+1:n,k)
U=triu(A);
y=zeros(n,1);
for k=1:n
y(k)=(b(k)-L(k,1:k-1)*y(1:k-1))/L(k,k);
end
x=zeros(n,1);
x(n)=y(n)/U(n,n);
for k=n-1:-1:1
x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);
end

function [x,L,U,D]=xmycholesky(A,b)
n=length(b);
U=zeros(n,n);
L=eye(n,n);
D=zeros(1,n);
U(1,1:n)=A(1,1:n);
L(2:n,1)=U(1,2:n)/U(1,1);
for k=2:n
U(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);
L(k+1:n,k)=U(k,k+1:n)/U(k,k);
end
D=diag(diag(U));
y=zeros(n,1);
y(1)=b(1);
for k=2:n
y(k)=b(k)-L(k,1:k-1)*y(1:k-1);
end
for k=1:n
Z(k)=y(k)/D(k,k);
end
U=L';
x=zeros(n,1);
x(n)=Z(n);
for k=n-1:-1:1
x(k)=Z(k)-U(k,k+1:n)*x(k+1:n);
end

function det=xmydet(A)
det=1;
n=length(A);
for k=1:n-1
maxvalue=max(abs(A(k:n,k:n)));
maxValue=max(maxvalue);
[ik,jk]=find(A==maxValue);
if ik~=k
A([k ik],:)=A([ik k],:);det=-det;
end
if jk~=k
A(:,[k jk])=A(:,[jk k]);det=-det;
end
A(k+1:n,k)=A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);
det=det*A(k,k);
end
det=det*A(n,n);

function [x]=xmychase(L,B,U,b)
n=length(B);
L=L';U=U';
for k=2:n
B(k)=B(k)-(L(k)/B(k-1))*U(k-1);
b(k)=b(k)-(L(k)/B(k-1))*b(k-1);
end
x(n)=b(n)/B(n);
for k=n-1:-1:1
x(k)=(b(k)-U(k)*x(k+1))/B(k);
end