数值分析中求解线性方程组的MATLAB程序(6种).doc

(完整word)数值分析中求解线性方程组的MATLAB程序(6种)
数值分析中求解线性方程组的MATLAB程序(6种）
1。回溯法（系数矩阵为上三角)
function X=uptrbk(A,B)
％求解方程组，首先化为上三角，再调用函数求解
[N，N］=size(A）；
X=zeros(N，1）；
C=zeros（1,N+1)；
Aug=[A B];
for p=1：N—1
［Y,j］=max（abs(Aug(p:N，p)));
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
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
D=Aug；
X=backsub(Aug（1:N，1：N），Aug（1:N，N+1));

2。系数矩阵为下三角
function x=matrix_down（A,b)
％求解系数矩阵是下三角的方程组
n=length（b);
x=zeros(n,1）;
x（1）=b（1）/A（1,1);
for k=2:1：n
x（k)=(b(k)-A（k,1：k-1)*x（1：k—1)）/A(k,k）;
end

3。普通系数矩阵（先化为上三角，在用回溯法）
function X=uptrbk(A,B)
％求解方程组，首先化为上三角,再调用函数求解
[N，N］=size(A）;
X=zeros(N,1）；
C=zeros（1,N+1)；
Aug=[A B];
for p=1:N-1
［Y,j］=max（abs(Aug(p:N,p)））；
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
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
D=Aug;
X=backsub（Aug(1:N，1：N)，Aug(1:N，N+1）)；

4.三角分解法
function ［X,L，U］=LU_matrix(A，B)
%A是非奇异矩阵
％AX=B化为LUX=B，L为下三角,U为上三角
%程序中并没有真正解出L和U，全部存放在A中
［N，N］=size(A);
X=zeros（N,1)；
Y=zeros(N，1）；
C=zeros（1,N）;
R=1:N；
for p=1:N-1
［max1，j]=max(abs(A（p：N,p)）);
C=A（p，:）；
A（p,：)=A（j+p-1，：);
A(j+p-1，:）=C；
d=R(p）；
R（p）=R（j+p-1）;
R(j+p—1)=d;
if A(p,p）==0
'A is singular.No unique solution'
break；
end
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
Y（1)=B（R（1)）；
for k=2：N
Y(k)=B(R(k））-A（k，1：k—1)＊Y（1:k—1);
end
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
L=tril（A，-1)+eye(N）
U=triu(A)

5.雅克比迭代法
function X=jacobi(A,B，P,delta，max1）;
％雅克比迭代求解方程组
N=length（B）；
for k=1：max1
for j=1:N
X（j)=(B(j)-A（j,［1：j-1，j+1:N］)*P([1：j—1，j+1：N］））/A（j,j)；
end
err=abs（norm(X’-P)）；
relerr=err/（norm(X)+eps)；
P=X’;
if （err〈delta）｜(relerr<delta)
break
end
end
X=X';
k

6.盖斯迭代法
function X=gseid(A，B，P，delta,max1);
％盖斯算法，求解赋初值的微分方程
N=length(B);
for k=1:max1
for j=1:N
if j==1
X（1)=（B(1)-A（1,2:N)*P（2：N))/A(1，1）;
elseif j==N
X(N）=(B(N）-A(N,1：N—1)*（X(1：N-1））’)/A(N，N）;
else
X（j)=(B（j)-A(j,1：j—1）*X（1:j-1)-A（j，j+1：N)*P(j+1:N））/A(j,j);
end
end
err=abs（norm（X'—P)）；
relerr=err/(norm（X)+eps);
P=X';
if (err<delta）|（relerr<delta)
break；
end
end
X=X';
k