帕斯卡中的赛德尔方法

如何解决帕斯卡中的赛德尔方法

我需要在 Pascal 中实现 Seidel 方法。我试过这段代码，但它给出了错误的答案。我不明白错误是什么。这是寻找根的过程：

CanDeactivate

变量'a'的过程

procedure Seidel(n: Integer; var x: vector; a: matrix; e: Real);
var k,i,j,z: integer;
s: Real;
begin
  for k := 1 to 100 do
    begin
      z := k;
      for i := 1 to n do
        begin
          s := a[i,n + 1];
          for j := 1 to n do s := s - a[i,j] * x[j];
          s := s / a[i,i];
          x[i] := x[i] + s;
          if abs(s) > e then z := 0
        end;
      if z <> 0 then Break;
    end;
end;

这是 StringGrid 的样子： "Корни Х" - "Roots X"

点击“Расчёт”（计算）按钮时，答案不同，反复点击后，出现“浮点溢出”错误。

解决方法

错误是

不使用评论
使用 2 个以上的单字母变量名
使用反模式：仅当您可以预测精确的迭代次数时，才应使用计数循环（for 循环）。 Break 确实/不应该属于您的标准曲目，我什至认为它是 spaghetti code 的变体。此规则很少有exception，但在这里您最好坚持使用条件循环（while … do 或 repeat … until）。
省略 begin … end 帧（用于分支和循环）开发过程中，当您的程序显然尚未完成时

公平地说，Seidel 方法可能令人困惑。另一方面，Pascal 提供了足够的语言能力，非常适合这样的任务。

实际上，我不得不自己编写该任务的程序，以便了解为什么您的 procedure 没有产生正确的结果。以下 program 使用了一些扩展的 Pascal (ISO 10206) 功能，例如模式和类型查询。为此，您需要一个符合 EP 的编译器，例如 GPC（GNU Pascal 编译器）。 AFAIK，Delphi 不支持这些功能，但解决任何缺陷应该很容易。

考虑到上述所有“错误”，您得出以下解决方案。

program seidel(output);

type
    naturalNumber = 1..maxInt value 1;

除非另有说明，否则以下所有 naturalNumber 值均使用 1 进行初始化。这是一个 EP 扩展。

    linearSystem(
            coefficientCount: naturalNumber;
            equationCount: naturalNumber
        ) = record
            coefficient: array[1..equationCount,1..coefficientCount] of real;
            result: array[1..coefficientCount] of real;
            solution: array[1..equationCount] of real;
        end;

当然，您可以根据主要使用场景对数据类型进行不同的构建。

{
    Approximates the solution of the passed linearSystem
    using the Gauss-Seidel method.
    
    system.solution should contain an estimate of the/a solution.
}
procedure approximateSolution(var system: linearSystem);
    { Returns `true` if any element along the main diagonal is zero. }
    { NB: There is a chance of false negatives. }
    function mainDiagonalNonZero: Boolean;
    var
        product: real value 1.0;
        n: naturalNumber;
    begin
        { Take the product of all elements along the main diagonal. }
        { If any element is zero,the entire product is zero. }
        for n := 1 to system.coefficientCount do
        begin
            product := product * system.coefficient[n,n];
        end;
        
        mainDiagonalNonZero := product <> 0.0;
    end;

此函数 mainDiagonalNonZero 用于提醒您可以在例程中“嵌套”例程。虽然它在下面只调用一次，但如果你像这样构造代码单元，它会稍微清理你的源代码。

type
    { This is more readable than using plain integer values. }
    relativeOrder = (previous,next);
var
    approximation: array[relativeOrder] of type of system.solution;

注意，approximation是在getNextApproximationResidual前面声明的，所以这个函数和approximateSolution的主块都可以访问相同的 向量。

    { Calculates the next approximation vector. }
    function getNextApproximationResidual: real;
    var
        { used for both,identifying the equation and a coefficient }
        n: naturalNumber;
        { used for identifying one term,i.e. coefficient × solution }
        term: 0..maxInt;
        { denotes a current error of this new/next approximation }
        residual: real;
        { denotes the largest error }
        residualMaximum: real value 0.0;
        { for simplicity,you could use `approximation[next,n]` instead }
        sum: real;
    begin
        for n := 1 to system.equationCount do
        begin
            sum := 0.0;
            
            for term := 1 to n - 1 do
            begin
                sum := sum + system.coefficient[n,term] * approximation[next,term];
            end;
            
            { term = n is skipped,because that's what we're calculating }
            
            for term := n + 1 to system.equationCount do
            begin
                sum := sum + system.coefficient[n,term] * approximation[previous,term];
            end;

很明显，您的实现不包含两个 for 循环。它不会迭代所有术语。

            sum := system.result[n] - sum;
            { everything times the reciprocal of coefficient[n,n] }
            approximation[next,n] := sum / system.coefficient[n,n];
            
            { finally,check for larger error }
            residual := abs(approximation[next,n] - approximation[previous,n]);
            if residual > residualMaximum then
            begin
                residualMaximum := residual;
            end;
        end;
        
        getNextApproximationResidual := residualMaximum;
    end;

我已将这个函数 getNextApproximationResidual 外包，以便我可以在下面的循环中编写更好的中止条件。

const
    { Perform at most this many approximations before giving up. }
    limit = 1337;
    { If the approximation improved less than this value,}
    { we consider the approximation satisfactory enough. }
    errorThreshold = 8 * epsReal;
var
    iteration: naturalNumber;
begin
    if system.coefficientCount <> system.equationCount then
    begin
        writeLn('Error: Gauss-Seidel method only works  ','on a _square_ system of linear equations.');
        halt;
    end;
    
    { Values in the main diagonal later appear as divisors,}
    { that means they must be non-zero. }
    if not mainDiagonalNonZero then
    begin
        writeLn('Error: supplied linear system contains ','at least one zero along main diagonal.');
        halt;
    end;

不要不要信任用户输入。在我们计算任何东西之前，请确保 system 满足一些基本要求。 halt（不带任何参数）是一个 EP 扩展。某些编译器的 halt 还接受 integer 参数以将错误条件传达给操作系统。

    { Take system.solution as a first approximation. }
    approximation[next] := system.solution;
    
    repeat
    begin
        iteration := iteration + 1;
        { approximation[next] is overwritten by `getNextApproximationError` }
        approximation[previous] := approximation[next];
    end
    until (getNextApproximationResidual < errorThreshold) or_else (iteration >= limit);

or_else 运算符是 EP 扩展名。它明确表示“懒惰/快捷评估”。这里没有必要，但我喜欢它。

    { Emit a warning if the previous loop terminated }
    { because of reaching the maximum number of iterations. }
    if iteration >= limit then
    begin
        writeLn('Note: Maximum number of iterations reached. ','Approximation may be significantly off,','or it does not converge.');
    end;
    
    { Finally copy back our best approximation. }
    system.solution := approximation[next];
end;

我将以下内容用于测试目的。 protected (EP) 对应于 Delphi 中的 const（我猜）。

{ Suitable for printing a small linear system. }
procedure print(protected system: linearSystem);
const
    totalWidth = 8;
    fractionWidth = 3;
    times = ' × ';
    plus = ' + ';
var
    equation,term: naturalNumber;
begin
    for equation := 1 to system.equationCount do
    begin
        write(system.coefficient[equation,1]:totalWidth:fractionWidth,times,system.solution[1]:totalWidth:fractionWidth);
        for term := 2 to system.coefficientCount do
        begin
            write(plus,system.coefficient[equation,term]:totalWidth:fractionWidth,system.solution[term]:totalWidth:fractionWidth);
        end;
        
        writeLn('⩰  ':8,system.result[equation]:totalWidth:fractionWidth);
    end;
end;

以下示例线性方程组采用 from Wikipedia，所以我“知道”了正确的结果：

{ === MAIN ============================================================= }
var
    example: linearSystem(2,2);
begin
    with example do
    begin
        { first equation }
        coefficient[1,1] :=  16.0;
        coefficient[1,2] :=   3.0;
        result[1] := 11.0;
        
        { second equation }
        coefficient[2,1] :=   7.0;
        coefficient[2,2] := -11.0;
        result[2] := 13.0;
        
        { used as an estimate }
        solution[1] := 1.0;
        solution[2] := 1.0;
    end;
    
    approximateSolution(example);
    
    print(example);
end.