非线形方程组的解法，100分(100分)

最好是拟牛顿法和代码
欢迎更多的好建议

难道就没人帮我解决吗？
如果谁先解决，100分全给。

type
TFuncGroup = procedure (x: array of Double
var y: array of Double);
TMatrix = array of array of Double;

// 用差商代替导数的牛顿法
// 求非线性方程组
// f(x[1], x[2], ..., x[n]) = 0 (i = 1, 2, ..., n) (1)
//的实根 (x[1]*, x[2]*, ..., x[n]*)。
// 给定 ||f|| (i = 1, 2, ..., n) 所允许的误差 eps > 0，以初始向量 X[0, 0]
//= (x[0, 1], x[0, 2], ..., x[0, n]) 出发进行迭代，如第 k 次迭代值
// X[k, 0] = (x[k, 1], x[k, 2], ..., x[k, n])
//满足 ||f(X[k, 0])|| < eps (i = 1, 2, ..., n)
//则取 X[k, 0] 作为(1)的近似解；否则，构造几个新点
// X[k, i] = (x[k, 1], ..., x[k, i - 1], x[k, i] + h[k], ..., x[k, n]) (i = 1, 2, ..., n)
//其中 h[k] 是第 k 次迭代中的步长。然后对每一 f 作满足
// w[k, i](X[k, j]) = f(X[k, j]) (j = 0, 1, ..., n
i = 1, 2, ..., n)
//的一次内插多项式
// w[k, i](x[1], x[2], ..., x[n])
// = f(X[k, 0]) + Sum((f(X[k, j]) - f(X[k, 0]))(x[j] - x[k, j]) / h[k]) (i = 1, 2, ..., n)
//把方程
// w[k, i](x[1], x[2], ..., x[n]) = 0 (i = 1, 2, ..., n) (2)
//的解作为(1)的解的第 k + 1 次逼近。
// 为解(2)，作仿射变换。令
// a[k] = 1 - Sum(Z[k, m]), x[j] = x[k + 1, j] = x[k, j] - Z[k, j] h[k] / a[k] (j = 1, 2, ..., n)
//则解(2)可化为解
// Sum(f(X[k, j]) Z[k, j]) = f(X[k, 0]) (i = 1, 2, ..., n) (3)
//在 k -> Inf, {h[k]} -> 0时，如果{X[k]}收敛，则收敛于(1)的解。
//
//n 方程组阶数
//h 步长参数
//w 每次迭代中用以乘 h 的因子
//ep 方程所允许的误差
//ps 方程左边所允许的最大绝对值
//x[1..n] 开始时存放初始向量，结束时存放近似解
//y[1..n] 存放方程组左面各方程在X上的值
//F(x, y) 对X计算方程组左边各方程的值并存放在y中
//
//返回值 -1: 如果(1)左边有一个方程的绝对值大于ep，改变初始向量X或其他参数后，重新计算
// -2: 如果方程组(3)是奇异，用返回的X重新计算
// -3: 如果(3)的根的和等于1，用返回的X重新计算
// m: 达到精度要求时的迭代次数
//
//Eg: x[1] ^ 2 + x[2] ^ 2 - 1 = 0
// 0.75 * x[1] ^ 3 + 0.9 - x[2] = 0
//
// X = {-7, -2}, ep = 1E-7, w = 0.1, ps = 1E4, h = 1
// Res: (-0.9817026, 0.1904204}

function NLS(n: Integer
h, w, ep, ps: Double
var x, y: array of Double
F: TFuncGroup): Integer;

implementation

function NLS(n: Integer
h, w, ep, ps: Double
var x, y: array of Double
F: TFuncGroup): Integer;
var
m, i, k: Integer;
r, alpha: Double;
b1, b2: Boolean;
a: TMatrix;

function Gauss(u: Integer
a: TMatrix
var y: array of Double): Integer;
var
te: Double;
i, j, k, m, n: Integer;
Err: Boolean;
begin
n := 0;
Result := 0;
repeat
Inc;
Err := True;
for k := n - 1 to u - 1 do
if Abs(a[k, n]) > 1E-8 then begin
Err := False;
Break;
end;
if Err then begin
Result := -2;
Exit;
end;
if k + 1 <> n then
for m := n - 1 to u do begin
te := a[n, m];
a[n, m] := a[k, m];
a[k, m] := te;
end;
for j := u downto n - 1 do
a[n - 1, j] := a[n - 1, j] / a[n - 1, n - 1];
for i := k + 1 to u - 1 do
for j := n to u do
a[i, j] := a[i, j] - a[i, n - 1] * a[n - 1, j];
until n + 1 = u;
for i := u - 1 downto 0 do begin
y := a[i, u] / a[i, i];
for k := i - 1 downto 0 do
a[k, u] := a[k, u] - a[k, i] * y;
end;
end;

begin
SetLength(a, n, n + 1);
m := 0;
repeat
b1 := True;
b2 := False;
F(x, y);
for i := 0 to n - 1 do begin
r := y;
a[i, n] := r;
r := Abs(r);
b1 := b1 and (r < ep);
b2 := b2 or (r > ps);
end;
if b1 then begin
Result := m;
Exit;
end;
if b2 then begin
Result := -1;
Exit;
end;
for i := 0 to n - 1 do begin
r := x;
x := r + h;
F(x, y);
for k := 0 to n - 1 do
a[k, i] := y[k];
x := r;
end;
Result := Gauss(n, a, y);
if Result < 0 then Exit;
alpha := 1;
for i := 0 to n - 1 do
alpha := alpha - y;
if alpha = 0 then begin //
Result := -3;
Exit;
end;
alpha := h / alpha;
for i := 0 to n - 1 do
x := x - y * alpha;
h := h * w;
Inc(m);
until False;
end;

啊，忘了例子。

procedure Func(x: array of Double
var y: array of Double);
begin
y[0] := Sqr(x[0]) + Sqr(x[1]) - 1;
y[1] := 0.75 * Sqr(x[0]) * x[0] + 0.9 - x[1];
end;

procedure TForm1.FormClick(Sender: TObject);
var
x, y: array [0..1] of Double;
r: Integer;
begin
x[0] := -7
x[1] := -2;
y[0] := 0
y[1] := 0;
r := NLS(2, 1, 0.1, 1E-7, 1E4, x, y, Func);
ShowMessage(Format('%d: (%s, %s)',
[r, FormatFloat('0.0000000', x[0]), FormatFloat('0.0000000', x[1])]));
end;

谢谢JohnsonGuo大侠，能留个email吗？以后会有问题请教。
我的email是weisansao@263.net

johnson_guo0215@sina.com