发布时间 : 星期日 文章王能超 计算方法——算法设计及MATLAB实现课后代码更新完毕开始阅读2b5c7cbbad02de80d5d84004
x=x1;
算例16：用对称超松弛(SSOR)迭代法和超松弛(SSOR)迭代法求解线性方程组：
?4x1?x2?1? ??x1?4x2?x3?4
??x?4x??33?2精度控制在10.
解：令A=[4,-1,0;-1,4,-1;0,-1,4]; b=[1,4,-3]’; x0=[0,0,0]’;
在命令窗口中输入：
[x1,k1]=SORmethod(A,b,x0,100,10^-5,1.2) [x2,k2]=SSORmethod(A,b,x0,100,10^-5,1.2)
第六章 线性方程组的直接法 6.1 追赶法
用追赶法解三对角线性方程组Ax?f. MATLAB文件：（threedia.m） function x=threedia(a,b,c,f)
%求解线性方程组Ax=f，其中A是三对角阵 %a是矩阵A的下对角线元素a(1)=0 %b是矩阵A的对角线元素
%c是矩阵A的上对角线元素c(N)=0 %f是方程组的右端向量 N=length(f);
x=zeros(1,N); y=zeros(1,N); d=zeros(1,N); u=zeros(1,N); %预处理 d(1)=b(1); for i=1:N-1
u(i)=c(i)/d(i);
d(i+1)=b(i+1)-a(i+1)*u(i); end
%追的过程 y(1)=f(1)/d(1); for i=2:N
y(i)=(f(i)-a(i)*y(i-1))/d(i); end
%赶的过程 x(N)=y(N); for i=N-1:-1:1
x(i)=y(i)-u(i)*x(i+1);
21
?5end
算例17：用追赶法求解方程组：
0??x1??6??2?10??13?20??x??1???2???? ??0?12?1??x3??0???????00?35???x4??1?解：令a=[0,-1,-1,-3]; b=[2,3,2,5]; c=[-1,-2,-1,0]; f=[6,1,0,1]’;
在命令窗口输入：
x=threedia(a,b,c,f)
6.2 Cholesky方法
用Cholesky分解法解对称方程组Ax?b. MATLAB文件：（文件名：Chol_decompose.m） function x=Chol_decompose(A,b) %用Cholesky分解求解线性方程组Ax=b %A是对称矩阵
%L是单位下三角矩阵 %D是对角阵
%对矩阵A进行三角分解:A=LDL’ N=length(A);
L=zeros(N,N); D=zeros(1,N); for i=1:N L(i,i)=1; end
D(1)=A(1,1); for i=2:N for j=1:i-1 if j==1
L(i,j)=A(i,j)/D(j); else
sum1=0; for k=1:j-1
sum1=sum1+L(i,k)*D(k)*L(j,k); end
L(i,j)=(A(i,j)-sum1)/D(j); end end sum2=0; for k=1:i-1
sum2=sum2+L(i,k)^2*D(k); end
22
D(i)=A(i,i)-sum2; end
%分别求解线性方程组Ly=b；L’x=y/D y=zeros(1,N); y(1)=b(1); for i=2:N sumi=0; for k=1:i-1
sumi=sumi+L(i,k)*y(k); end
y(i)=b(i)-sumi; end x=zeros(1,N); x(N)=y(N)/D(N); for i=N-1:-1:1 sumi=0; for k=i+1:N
sumi=sumi+L(k,i)*x(k); end
x(i)=y(i)/D(i)-sumi; end
算例18：用Cholesky方法求解线性方程组：
?4?24??x1??8.7??????? ?21710x2?13.7 ????????4109????x3?????0.7??解：令A=[4,-2,4;-2,17,10;4,10,9]; b=[8.7,13.7,-0.7];
在命令窗口输入：
x=Chol_decompose(A,b)
6.3 矩阵分解方法
基于Gauss消去法的LU分解求解线性方程组Ax?b. MATLAB文件：（文件名：lu_decompose.m） function x=lu_decompose(A,b)
%基于Gauss消去法的LU分解求解线性方程组Ax=b %A表示系数矩阵
%b表示方程组右边的向量 n=length(b);
L=eye(n); U=zeros(n,n); x=zeros(n,1); y=zeros(n,1);
%对矩阵A进行LU分解
23
for i=1:n
U(1,i)=A(1,i); if i==1
L(i,1)=1; else
L(i,1)=A(i,1)/U(1,1); end end for i=2:n for j=i:n sum=0; for k=1:i-1
sum=sum+L(i,k)*U(k,j); end
U(i,j)=A(i,j)-sum; if j~=n sum=0; for k=1:i-1
sum=sum+L(j+1,k)*U(k,i); end
L(j+1,i)=(A(j+1,i)-sum)/U(i,i); end end end
%解方程组Ly=b y(1)=b(1); for k=2:n sum=0; for j=1:k-1
sum=sum+L(k,j)*y(j); end
y(k)=b(k)-sum; end
%解方程组Ux=y x(n)=y(n)/U(n,n); for k=n-1:-1:1 sum=0; for j=k+1:n
sum=sum+U(k,j)*x(j); end
x(k)=(y(k)-sum)/U(k,k); end
24