容差系统的Jacobi方法

Question

晚上所有，我以前有过类似的问题，这种类型的代码，但我不能弄清楚这一点。我正在尝试运行Jacobi代码，初始近似为0向量，容差矩阵范数(X^n - x^(n-1)) < 1e^-2

restart;

jacobi:=proc(A::Matrix,b::Vector,x0::Vector,m::integer)
  local n,i,j,S,x::Vector,x1::Vector; 
  x1:=Vector(x0); x := Vector(m);
ee:= LinearAlgebra[VectorNorm](x-x1);
  for n from 1 to 50 while ee < 1e-2 do
    for i from 1 to m do
       S:=0.0; 
       for j from 1 to m do
           if j<>i then S:=S+A[i,j]*x1[j]; end if;    
       end do;
       x[i]:=(-S+b[i])/A[i,i]; 
    end do;
    for j from 1 to m do 
       x1[j]:=x[j]; 
    end do;
    print(n,x);
  end do;
  return x
end proc;

    A:=Matrix([[12.4,0.3,-2.2,0.2,3.6],[1.2,7.1,-1.7,-1.6,0.9],[1.3,3.1,10.8,2.2,0.7],[-1.4,0.8,1.1,-7.7,-1.8],[2.8,6.4,0.0,-1.2,-16.6]]);

b:=Vector([-3.5,19.58,-24.56,8.16,-9.28]);


x0:=Vector([0,0,0,0,0]);

jacobi(A,b,x0,5);

#Error, (in jacobi) Vector index out of range

Error, (in LinearAlgebra:-Norm) expects its 1st argument, A, to be of type {Matrix, Vector}, but received x-x1

Answer 1

我找到了一个解决方案！

restart;
jacobi:=proc(A::Matrix,b::Vector,x0::Vector,m::integer,N::integer,tol::float)
 local n,i,j,S,ee,x::Vector,x1::Vector;
   x1:=Vector(x0); x := Vector(m); n:=1;
   do
      for i from 1 to m do
         S:=0.0; 
         for j from 1 to i-1 do
             S:=S+A[i,j]*x1[j]
         end do;
         for j from i+1 to m do
             S:=S+A[i,j]*x1[j]
         end do;
         x[i]:=(-S+b[i])/A[i,i];  
      end do;
      ee:= LinearAlgebra[VectorNorm](x-x1);
      if(ee<tol) then  
         printf("The solution after %d iterations is:",n);
         print(x); 
         break; 
      end if;
      for j from 1 to m do 
         x1[j]:=x[j]; 
      end do;
      n:=n+1;
      if (n>N) then error("No convergence after",N,"iterations") end if;  
   end do;
   return x;
 end proc;
 jacobi := proc(A::Matrix, b::Vector, x0::Vector, m::integer, 
   N::integer, tol::float)
    ... 
   end;



    A:=Matrix([[12.4,0.3,-2.2,0.2,3.6],[1.2,7.1,-1.7,-1.6,0.9],[1.3,3.1,10.8,2.2,0.7],[-1.4,0.8,1.1,-7.7,-1.8],[2.8,6.4,0.0,-1.2,-16.6]]);

b:=Vector([-3.5,19.58,-24.56,8.16,-9.28]);

x0:=Vector([0,0,0,0,0]);
                              
jacobi(A,b,x0,5,100,0.01);

经过8次迭代后的解决方案是：