disp('原方程为AX=b:') A b
disp('------------------------') n=length(b);
for k=1:n-1
[main,index]=max(abs(A(k:n,k))); index=index+k-1; if abs(main) disp('列元素太小!!'); break; elseif index>k temp=A(k,:); A(k,:)=A(index,:); A(index,:)=temp; end for i=k+1:n m(i,k)=A(i,k)/A(k,k); A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n); b(i)=b(i)-m(i,k)*b(k); end end disp('消元后所得到的三角阵是') A b(n)=b(n)/A(n,n); for i=n-1:-1:1 b(i)=(b(i)-sum(A(i,i+1:n).*b(i+1:n)))/A(i,i); end clear x; x=b; disp('AX=b的解x是') x 5 (二)实验题目:道立特三角分解法解方程组 一、程序功能 编程实现三角分解法解方程组。 二、实验算例选择 x1+2x2+x3=2 -2x1-2x2-x3=-3 2x1-3x2-2x3=-1 三、算法 步骤1 : 计算L下三角矩阵; 步骤2: 计算U; 步骤3: Ly=b,Ux=y求出x,y。 四、重要标识符说明 RA代表A的秩 RB代表B的秩 五、程序运行实例 六、源程序 A=[1 2 1;-2 -2 -1;2 -3 -2]; b=[2 ;-3 ;-1]; B=[A b]; n=length(b); RA=rank(A); RB=rank(B); 6 zhica=RB-RA; if zhica>0 disp('请注意:因为RA~=RB,所以此方程组无解.') return; end if RA==RB if RA==n disp('请注意:因为RA=RB=n,所以此方程组有唯一解.') X=zeros(n,1); for p= 1:n-1 for k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1); end end b=B(1:n,n+1); A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end disp('AX=b的解x是') X else disp('请注意:因为RA=RB 七、个人实验总结 高斯列主元消元法就是在高斯消元法的基础上选取每列最大值最为对角线上的元素,以减小误差;而三角分解法需求出上三角和下三角然后解出X,Y并不用选主元。但前者数值较为稳定。 7 (三)实验题目:雅可比解方程组 一、程序功能 编程实现雅可比法解方程组 二、实验算例选择 A=[-4,1,1,1;1,-4,1,1,1;1,1,-4,1;1,1,1,-4] B=[1;1;1;1] ; 三、算法 步骤1:求出矩阵的范数,判断迭代公式是否收敛。 步骤2: 代入雅克比的矩阵迭代公式x0迭代,得到解 四、重要标识符说明 A为方程组的系数矩阵 b为方程组的常数向量 X0为迭代初值 Eps为误差限 T为最大迭代次数,超过则认为无解。 五、程序运行实例 六、源程序 function jacobi clc; clear; A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 0 0 -1;-1 0 0 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; 8 B=[0 ;5 ;0 ;6 ;-2 ;6]; Err_user=0.01; N=50; [m,n]=size(A); X=zeros(n,1); k=1; while k<=N Xk=X; for i=1:n for j=1:n if i~=j AX(j)=A(i,j)*Xk(j); end end Sum_AX=sum(AX); AX=0; X(i)=(B(i)-Sum_AX)/A(i,i); end E=max(abs(Xk-X)); if E (三)实验题目:高斯—赛德尔解方程组 一、程序功能 编程实现雅可比法解方程组 二、实验算例选择 A=[-4,1,1,1;1,-4,1,1,1;1,1,-4,1;1,1,1,-4] B=[1;1;1;1] 三、算法 步骤1 :求出矩阵的范数,判断迭代公式是否收敛。 9