线性方程组数值方法

线性方程组数值方法
Ax ，必须满足系数矩阵的对角元素非零的条件。

1、回代：用回代法求解上三角线性方程组b
function x=backsub(A,b)
%Input -A is an nxn upper-triangular nonsingular matrix
% -b is an nx1 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=[1 2 1 4;0 -4 2 -5;0 0 -5 -7.5;0 0 0 -9];b=[13;2;-35;-18];
x=backsub(A,b)
2、高斯消去法(不选主元素)
function x=gauss(A,b)
%Input -A is an nxn nonsingular matrix
% -b is an nx1 matrix
%Output-x is an nx1 matrix containing the solution to Ax=b
%Initialize x
[n,n]=size(A);x=zeros(n,1);
%Form the augmented matrix:Aug=[A|b]
Aug=[A b];
for k=1:n-1
%Elimination process for column k
for i=k+1:n
m=Aug(i,k)/Aug(k,k);
Aug(i,k:n+1)=Aug(i,k:n+1)-m*Aug(k,k:n+1);
end
end
%Output Elimination matrix
Aug
%Back Substitution
A=Aug(1:n,1:n);b=Aug(1:n,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=[4 -9 2;2 -4 6;-1 1 -3];b=[5;3;-4]; x=gauss(A,b)
3、列主元消去法
function x=gauss_column(A,b)
%Input -A is an nxn nonsingular matrix
% -b is an nx1 matrix
%Output-x is an nx1 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 k=1:n-1
%Partial pivoting for column k
[y,j]=max(abs(Aug(k:n,k)));
%Interchange row k and j
C=Aug(k,:);Aug(k,:)=Aug(j+k-1,:);Aug(j+k-1,:)=C;
if Aug(k,k)==0
'A was singular,No unique solution'
break
end
%Elimination process for column k
for i=k+1:n
m=Aug(i,k)/Aug(k,k);
Aug(i,k:n+1)=Aug(i,k:n+1)-m*Aug(k,k:n+1);
end
end
%Back Substitution on [U|y] using backsub
x=backsub(Aug(1:n,1:n),Aug(1:n,n+1));
例：A=[1 2 1 4;2 0 4 3;4 2 2 1;-3 1 3 2];b=[13;28;20;6];
x=gauss_column(A,b)
4、LU分解法(不选主元素)
function x=lufact(A,b)
%Input -A is an nxn matrix
% -b is an nx1 matrix
%Output-x is an nx1 matrix containing the solution to Ax=b
%Initilize x,y,
[n,n]=size(A);x=zeros(n,1);y=zeros(n,1);
for k=1:n-1
%Calculate multiplier and place in subdiagonal portion of A for i=k+1:n
mult=A(i,k)/A(k,k);
A(i,k)=mult;
A(i,k+1:n)=A(i,k+1:n)-mult*A(k,k+1:n);
end
end
%Output LU factorization compact matrix
A
%solve for y
y(1)=b(1);
for i=2:n
y(i)=b(i)-A(i,1:i-1)*y(1:i-1);
end
%Solve for x
x(n)=y(n)/A(n,n);
for i=n-1:-1:1
x(i)=(y(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);
end
例：A=[4 -9 2;2 -4 6;-1 1 -3];b=[5;3;-4];x=lufact(A,b)
5、列主元LU分解法（PA＝LU）A是非奇异矩阵
function x=lu_column(A,b)
%Input -A is an nxn matrix
% -b is an nx1 matrix
%Output-x is an nx1 matrix containing the solution to Ax=b
%Initilize x,y,the temporary storage matrix C,and the row
%permutation information matrix R
[n,n]=size(A);
x=zeros(n,1);y=zeros(n,1);C=zeros(1,n);R=1:n;
for k=1:n-1
%Find the pivot row for column k
[maxl,j]=max(abs(A(k:n,k)));
%Interchange row k and j
C=A(k,:);A(k,:)=A(j+k-1,:);A(j+k-1,:)=C;
d=R(k);R(k)=R(j+k-1);R(j+k-1)=d;
if A(k,k)==0
'A is singular,No unique solution'
break
end
%Calculate multiplier and place in subdiagonal portion of A for i=k+1:n
mult=A(i,k)/A(k,k);
A(i,k)=mult;
A(i,k+1:n)=A(i,k+1:n)-mult*A(k,k+1:n);
end
end
%Output LU factorization compact matrix
A
%solve for y
y(1)=b(R(1));
for i=2:n
y(i)=b(R(i))-A(i,1:i-1)*y(1:i-1);
end
%Solve for x
x(n)=y(n)/A(n,n);
for i=n-1:-1:1
x(i)=(y(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);
end
例：A=[1 2 6;4 8 -1;-2 3 -5];b=[1;2;3];
x=lu_column(A,b)
6、Jacobi迭代:程序可用的充分条件是A具有严格对角优势。

function [x,k,err,relerr]=jacobi(A,b,x0,delta,maxl)
%Input -A is an nxn nonsingular matrix
% -b is an nx1 matrix
% -x0 is an nx1 matrix;the initial guess
% -delta is the tolerance for x0
% -maxl is the maximum number of iterations
%Output-x is an nx1 matrix;the jacobi approximation to the
% solution of Ax=b
n=length(b);x=zeros(n,1);
for k=1:maxl
for i=1:n
x(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i); end
err=abs(norm(x-x0));
relerr=err/(norm(x)+eps);
x0=x;
if(err<delta)|(relerr<delta)
break
end
end
例：A=[4 -1 1;4 -8 1;-2 1 5];b=[7;-21;15];
[x,k,err,relerr]=jacobi(A,b,[1;2;2],0.001,50)
7、Gauss-Seidel迭代:程序可用的充分条件是A具有严格对角优势。

function [x,k,err,relerr]=gauss_seidel(A,b,x0,delta,maxl)
%Input -A is an nxn nonsingular matrix
% -b is an nx1 matrix
% -x0 is an nx1 matrix;the initial guess
% -delta is the tolerance for x0
% -maxl is the maximum number of iterations
%Output-x is an nx1 matrix;the gauss_seidel approximation to the
% solution of Ax=b
n=length(b);x=zeros(n,1);
for k=1:maxl
for i=1:n
if i==1
x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);
elseif i==n
x(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);
else
%x contains the kth approximations and x0 the (k-1)st
x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i); end
end
err=abs(norm(x-x0));
relerr=err/(norm(x)+eps);
x0=x;
if(err<delta)|(relerr<delta)
break
end
end
例：A=[4 -1 1;4 -8 1;-2 1 5];b=[7;-21;15];
[x,k,err,relerr]=gauss_seidel(A,b,[1;2;2],0.001,50)
例1。