Метод Ньютона

[lordsyrius]

Здравствуйте. Написал программу, в которой есть нахождение корня методом Ньютона, но в чем-то допустил ошибку, т.к. ответ не правильный... помогите, пожалуйста.

Function f1(x0, y1, y2,y3,x12,x21,x13,x31:Real):real;
begin
f1:=1-y1-3*x0+x0*y1+3*x0*x0-x0*x0*x0+x12*x21-x12*x21*x0-x12*y2+x13*y3-x13*x31+x13*x31*x0;
end;

Function f2(x0,y1,x12,x21,x13,x31:Real):real;
begin
f2:=-3+y1+6*x0-3*x0*x0-x12*x21+x13*x31;
end;

procedure TForm1.btn1Click(Sender: TObject);
label
1, 2;
  var
x12,x13,x21,x23,x31,x32: real;
b1, b2, b3, b: real;
m1, m2, m3: real;
k1, k2, k3, L, s, p: real;
a1, a2, a3: real;
y1,y2,y3: real;
z,d, eps,x0, x,xn1,y0: real;
n:Integer;
begin
x12:=StrToFloat(edt2.Text);
x13:=StrToFloat(edt3.Text);
x21:=StrToFloat(edt4.Text);
x23:=StrToFloat(edt6.Text);
x31:=StrToFloat(edt7.Text);
x32:=StrToFloat(edt8.Text);

b1:= Exp((1/3)*ln(x12*x13));
b2:= Exp((1/3)*ln(x21*x23));
b3:= Exp((1/3)*ln(x31*x32));
b:=b1+b2+b3;

a1:=b1/b;
a2:=b2/b;
a3:=b3/b;

Label8.Caption:= 'Вектор приоритетов:';
Label5.Caption:= FloatToStrF(a1,ffFixed,5,5);
Label6.Caption:= FloatToStrF(a2,ffFixed,5,5);
Label7.Caption:= FloatToStrF(a3,ffFixed,5,5);
img3.Visible:=True;
img4.Visible:=True;

m1:=1*a1+x12*a2+x13*a3;
m2:=x21*a1+1*a2+x23*a3;
m3:=x31*a1+x32*a2+1*a3;

k1:=m1/a1;
k2:=m2/a2;
k3:=m3/a3;

L:=(k1+k2+k3)/3;
s:=(L-3)/2;
p:=s/0.58;
lbl5.Caption:= ' лямда max = ' + FloatToStrF(L,ffFixed,5,5);
Label1.Caption:= ' Индекс согласованности = ' + FloatToStrF(s,ffFixed,5,5);
Label2.Caption:= ' Отношение согласованности = ' + FloatToStrF(p,ffFixed,5,5);

y1:=x23*x32;
y2:=x23*x31;
y3:=x21*x32;

Label4.Caption:= 'Характеристическое уравнение: ';
Label3.Caption:= '(1-L) * [(1-L)^2 - '+FloatToStrF(y1,ffFixed,4,2)+'] + '+ FloatToStrF(x12,ffFixed,4,2)+' * ['+FloatToStrF(x21,ffFixed,4,2)+' * (1-L) - '+FloatToStrF(y2,ffFixed,4,2)+'] + '+FloatToStrF(x13,ffFixed,4,2)+' * ['+FloatToStrF(y3,ffFixed,4,2)+' - '+FloatToStrF(x31,ffFixed,4,2)+' * (1-L)] = 0';

////////////////////////////Метод Ньютона//////////////////////////

x:=5;
eps:=0.01;
n:=0;

1: x0:=x;
if f2(x0,y1,x12,x21,x13,x31)<>0 then
begin
x:=x0-f1(x0, y1, y2,y3,x12,x21,x13,x31)/f2(x0,y1,x12,x21,x13,x31);
2: n:=n+1;
while Abs(x0-x)>eps do  goto 1;
end
else goto 2;
Label9.caption:=FloatToStrF(x,ffFixed,4,4)+'Число итераций ' + IntToStr(n);

В программе пользователь вводит симметричную матрицу и жмет на кнопку, после чего выводятся лямда мах, индекс согласованности, отношение согласованности, характеристическое уравнение, и должен быть найден по этому характеристическому уравнению методом Ньютона корень, который должен, по идеи, равняться лямде мах, которая найдена ранее.

[lordsyrius]

Неужели никто не сталкивался с методом Ньютона!?

Для лучшего понимания, написал комментарии к коду:

Function f1(x0, y1, y2,y3,x12,x21,x13,x31:Real):real; //наша функция, выведенная из 
                                                      //характеристического уравнения ниже
begin
f1:=1-y1-3*x0+x0*y1+3*x0*x0-x0*x0*x0+x12*x21-x12*x21*x0-x12*y2+x13*y3-x13*x31+x13*x31*x0;
end;

Function f2(x0,y1,x12,x21,x13,x31:Real):real;        // производная от нашей функции
begin
f2:=-3+y1+6*x0-3*x0*x0-x12*x21+x13*x31;
end;

procedure TForm1.btn1Click(Sender: TObject);    
label
1, 2;
  var
x12,x13,x21,x23,x31,x32: real;
b1, b2, b3, b: real;
m1, m2, m3: real;
k1, k2, k3, L, s, p: real;
a1, a2, a3: real;
y1,y2,y3: real;
z,d, eps,x0, x,xn1,y0: real;
n:Integer;

{расчет лямды мах (наш корень) методом операций с матрицами + расчет  остальных показатей (в том числе y1, y2, y3 для нашей функции)}
begin
x12:=StrToFloat(edt2.Text);  // элементы матрицы, которые задает пользователь (по диагонали единицы)
x13:=StrToFloat(edt3.Text);
x21:=StrToFloat(edt4.Text);
x23:=StrToFloat(edt6.Text);
x31:=StrToFloat(edt7.Text);
x32:=StrToFloat(edt8.Text);

b1:= Exp((1/3)*ln(x12*x13));
b2:= Exp((1/3)*ln(x21*x23));
b3:= Exp((1/3)*ln(x31*x32));
b:=b1+b2+b3;

a1:=b1/b;
a2:=b2/b;
a3:=b3/b;

Label8.Caption:= 'Вектор приоритетов:';
Label5.Caption:= FloatToStrF(a1,ffFixed,5,5);
Label6.Caption:= FloatToStrF(a2,ffFixed,5,5);
Label7.Caption:= FloatToStrF(a3,ffFixed,5,5);
img3.Visible:=True;
img4.Visible:=True;

m1:=1*a1+x12*a2+x13*a3;
m2:=x21*a1+1*a2+x23*a3;
m3:=x31*a1+x32*a2+1*a3;

k1:=m1/a1;
k2:=m2/a2;
k3:=m3/a3;

L:=(k1+k2+k3)/3;
s:=(L-3)/2;
p:=s/0.58;

lbl5.Caption:= ' лямда max = ' + FloatToStrF(L,ffFixed,5,5);
Label1.Caption:= ' Индекс согласованности = ' + FloatToStrF(s,ffFixed,5,5);
Label2.Caption:= ' Отношение согласованности = ' + FloatToStrF(p,ffFixed,5,5);

y1:=x23*x32;
y2:=x23*x31;
y3:=x21*x32;

Label4.Caption:= 'Характеристическое уравнение: ';
Label3.Caption:= '(1-L) * [(1-L)^2 - '+FloatToStrF(y1,ffFixed,4,2)+'] + '+ FloatToStrF(x12,ffFixed,4,2)+' * ['+FloatToStrF(x21,ffFixed,4,2)+' * (1-L) - '+FloatToStrF(y2,ffFixed,4,2)+'] + '+FloatToStrF(x13,ffFixed,4,2)+' * ['+FloatToStrF(y3,ffFixed,4,2)+' - '+FloatToStrF(x31,ffFixed,4,2)+' * (1-L)] = 0'; ////////наша функция

////////////////////////////Метод Ньютона//////////////////////////

x:=5;  {как я понял из теории, это начальная точка, которую мы сами задаем и от неё уже движимся влево}
eps:=0.01;
n:=0;

1: x0:=x;
if f2(x0,y1,x12,x21,x13,x31)<>0 then
begin
x:=x0-f1(x0, y1, y2,y3,x12,x21,x13,x31)/f2(x0,y1,x12,x21,x13,x31);
2: n:=n+1;
while Abs(x0-x)>eps do  goto 1;
end
else goto 2;
Label9.caption:=FloatToStrF(x,ffFixed,4,4)+'Число итераций ' + IntToStr(n);

[Stilet]

А если пройтись по ссылкам внизу страницы?

[NiCola999]

while Abs(x0-x)>eps do  goto 1;

вы забыли еще одно условие

while Abs(x0-x)>eps And Abs(f1(x, y1, y2,y3,x12,x21,x13,x31)) > eps

кстати производную можно найти так:

F(x)' = (f(x + h) - f(x - h)) / (2 * h)
h примерно 0.0000001

[lordsyrius]

Спасибо, вроде правильно всё сделал...

[IvaniuS]

lordsyrius на будущее разберитесь со своими руками и заставите офформлять код получше, практически не читабельный.