错误是
- 不使用 cmets
- 使用超过 2 个单字母变量名
- 使用反模式：仅当您可以预测准确迭代次数时，才应使用计数循环（
for 循环）。 Break 确实/不应该属于您的标准曲目，我什至认为它是spaghetti code 的变体。这条规则很少有exceptions，但在这里你最好坚持使用条件循环（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 的主块都可以访问 same 向量。
{ 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,
times,
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.