王能超 计算方法——算法设计及MATLAB实现课后代码

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