Pascal 中的赛德尔法
Seidel method in Pascal
我需要在 Pascal 中实现 Seidel 方法。我试过这段代码，但它给出了错误的答案。我不明白错误是什么。这就是寻找根的过程：
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;
变量 'a'
的过程
procedure ReadA;
var i, j: integer;
begin
for i := 1 to m do
for j := 1 to m + 1 do
a[i, j] := StrToFloat(Form1.StringGrid1.Cells[j, i])
end;
这是 StringGrid 的样子：
"Корни Х" - "Roots X"
点击“Расчёт”（计算）按钮时，答案不一样，反复点击后，出现“浮点数溢出”错误
错误是
- 不使用注释
- 使用超过 2 个单字母变量名
- 使用反模式：仅当您可以预测准确[=90=时才应使用计数循环（
for循环） ] 迭代次数。 Break does/should 不属于你的标准曲目，我什至认为它是这个规则的 spaghetti code. There are very few exceptions 的变体，但在这里你最好坚持使用条件循环（while … do 或 repeat … until).
- 省略
begin … end 框架（用于分支和循环）在开发期间，当你的程序显然还没有完成时
公平地说，Seidel 方法 可能 令人困惑。另一方面，如果具备足够的语言能力，Pascal 非常适合此类任务。
我实际上必须自己编写该任务才能理解为什么您的 procedure 没有产生正确的结果。以下 program 使用了一些 Extended Pascal (ISO 10206) 功能，例如模式和类型查询。您将需要一个符合 EP 的编译器，例如 GPC（GNU Pascal 编译器）。据我所知，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的主块都可以访问the 相同 个向量。
{ 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 参数来将错误条件传达给 OS.
{ 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 扩展。它明确表示“lazy/short-cut 评估”。这里没有必要，但我还是喜欢。
{ 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.
我需要在 Pascal 中实现 Seidel 方法。我试过这段代码，但它给出了错误的答案。我不明白错误是什么。这就是寻找根的过程：
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;
变量 'a'
的过程procedure ReadA;
var i, j: integer;
begin
for i := 1 to m do
for j := 1 to m + 1 do
a[i, j] := StrToFloat(Form1.StringGrid1.Cells[j, i])
end;
这是 StringGrid 的样子： "Корни Х" - "Roots X"
点击“Расчёт”（计算）按钮时，答案不一样，反复点击后，出现“浮点数溢出”错误
错误是
- 不使用注释
- 使用超过 2 个单字母变量名
- 使用反模式：仅当您可以预测准确[=90=时才应使用计数循环（
for循环） ] 迭代次数。Breakdoes/should 不属于你的标准曲目，我什至认为它是这个规则的 spaghetti code. There are very few exceptions 的变体，但在这里你最好坚持使用条件循环（while … do或repeat … until). - 省略
begin … end框架（用于分支和循环）在开发期间，当你的程序显然还没有完成时
公平地说，Seidel 方法 可能 令人困惑。另一方面，如果具备足够的语言能力，Pascal 非常适合此类任务。
我实际上必须自己编写该任务才能理解为什么您的 procedure 没有产生正确的结果。以下 program 使用了一些 Extended Pascal (ISO 10206) 功能，例如模式和类型查询。您将需要一个符合 EP 的编译器，例如 GPC（GNU Pascal 编译器）。据我所知，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的主块都可以访问the 相同 个向量。
{ 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 参数来将错误条件传达给 OS.
{ 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 扩展。它明确表示“lazy/short-cut 评估”。这里没有必要，但我还是喜欢。
{ 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.