Міністерство освіти і науки України
Львівський національний університет ім. Івана Франка
Факультет прикладної математики та інформатики
Звіт
на тему:
“Метод Стеффенсена”
Розглядаючи метод Ньютона видно, що на кожній ітерації потрібно обраховувати матрицю других похідних. Звідси зрозуміло, що в тих випадках, коли обчислення матриці других похідних потребує значного об'єму обчислень, трудомісткість кожної ітерації методу Ньютона може бути обтяжуючою. Тому виникає питання, про можливість побудови методів мінімізації, які б по швидкості збігання не були гіршими від методу Ньютона, але для своєї реалізації не вимагали обчислення матриці других похідних. Одним з таких методів є метод Стеффенсена, який розроблявся для розв'язування нелінійних рівнянь. Використовуючи цей метод для розв'язування системи рівнянь , отримаємо наступний ітераційний метод розв'язування задачі мінімізації
(1)
Якщо наближення вже відомо, то наступне наближення визначається так:
(2)
де — числовий параметр, – матриця розділених різниць перших похідних, яка визначається за правилом
(3)
тут – елемент і - го рядка j - го стовпця матриці , a , позначають відповідно перші і другі похідні по змінних ui, uj функції J(u) i, j – (1,...,n); припускається, що при . Тоді з (3) випливає, що
, .
Це означає, що при , метод (2) перетворюється в метод Ньютона.
Як видно з (2), на кожному кроці методу Стеффенсена потрібно розв'язувати систему лінійних алгебраїчних рівнянь
(4)
Тут мається на увазі, що матриця не вироджена. Якщо при всіх J=1 … n, то згідно формули (3) для визначення елементів матриці достатньо знати перші похідні в точках u=uk і . Якщо ж при деякому j то в силу (3) для визначення j-ro стовпця матриці прийдеться обраховувати другі похідні в точці . Якщо при деякому k(0 виявилось J'( uk )=0, то процес (2) чи (4) припиняється: в цьому випадку uk– стаціонарна точка, і для вияснення, чи буде uk(U*, при необхідності потрібно провести додатковий експеримент. Тому будемо вважати, що . А тоді для визначення матриці потрібно буде обраховувати не всі другі похідні і в цьому випадку метод (2) має перевагу перед методом Ньютона. З іншої сторони, виявляється, метод (2) в класі сильно випуклих функції збігається з такою ж швидкістю, як і метод Ньютона.
Розглянемо реалізацію методу на основі функції
x*=(5,1,2)
x0=(0.5,0.5,0.5);
unit Unit1;
interface
uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, Grids, ExtCtrls, ComCtrls, StdCtrls, Buttons;
type
TForm1 = class(TForm)
StatusBar1: TStatusBar;
Panel1: TPanel;
StringGrid1: TStringGrid;
Panel2: TPanel;
SpeedButton1: TSpeedButton;
SpeedButton2: TSpeedButton;
procedure FormCreate(Sender: TObject);
Function f(x: Array of Double): Double;
Function fp1(x: Array of Double; c: Integer): Double;
procedure SpeedButton1Click(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;
Const
n=3;
var
Form1: TForm1;
implementation
{$R *.dfm}
procedure TForm1.FormCreate(Sender: TObject);
begin
StringGrid1.Cells[0, 0] := 'K';
StringGrid1.Cells[1, 0] := 'X1';
StringGrid1.Cells[2, 0] := 'X2';
StringGrid1.Cells[3, 0] := 'X3';
StringGrid1.Cells[5, 0] := 'f(x)';
end;
Function TForm1.f(x: Array of Double): Double;
begin
Result :=sqr(x[0]-5)*(x[0]-5)+sqr(x[1]-1)*sqr(x[1]-1)+sqr(x[2]-2)*sqr(x[2]-2);
end;
Function TForm1.fp1(x: Array of Double; c: Integer): Double;
begin
Case c of
0: Result := 3*sqr(x[0]-5);
1: Result := 4*sqr(x[1]-1)*(x[1]-1);
2: Result := 4*sqr(x[2]-2)*(x[2]-2);
end;
end;
procedure TForm1.SpeedButton1Click(Sender: TObject);
var
x, x0, xh, Vektor1, Vektor2: array of Double;
M: array[0..n-1, 0..n-1] of Double;
i, j, k, iter: Integer;
eps,s1,s2,s3: Double;
begin
SetLength(x,n);
SetLength(x0,n);
SetLength(Vektor1, n);
SetLength(Vektor2, n);
SetLength(xh, n);
x[0] := 0;
x[1] := 0;
x[2] := 0;
eps := 0.01;
iter := 0;
While (abs(f(x)) > eps)and(iter < 1000) do
begin
inc(iter);
x0[0] := x[0];
x0[1] := x[1];
x0[2] := x[2];
For i := 0 to n-1 do
xh[i] := x[i] - fp1(x,i);
For k := 0 to n-1 do
begin
For i := 0 to n-1 do
begin
Vektor1[i] := x[i];
Vektor2[i] := x[i];
end;
For i := 0 to k do
Vektor1[i] := xh[i];
For i := 0 to k - 1 do
Vektor2[i] := xh[i];
For i := 0 to n-1 do
if (xh[k] - x[k]) <> 0 then
M[i, k] := (fp1(Vektor1, i) - fp1(Vektor2,i))/(xh[k] - x[k])
else
M[i, k] := 0;
end;//for k
For i := 0 to n-1 do
begin
if M[i,0] <> 0 then
s1 := (1/M[i,0])*fp1(x,0)
else
s1 := 0;
if M[i,1] <> 0 then
s2 := (1/M[i,1])*fp1(x,1)
else
s2 := 0;
if M[i,2] <> 0 then
s3 := (1/M[i,2])*fp1(x,2)
else
s3 := 0;
x[i] := x[i] - (s1+s2+s3);
StringGrid1.Cells[i+1, StringGrid1.RowCount - 1] := FloatToStr(x[i]);
end;
StringGrid1.Cells[5, StringGrid1.RowCount - 1] := FloatToStr(f(x));
StringGrid1.Cells[0, StringGrid1.RowCount - 1] := FloatToStr(iter);
StringGrid1.RowCount := StringGrid1.RowCount + 1;
end;
end;
end.
Висновок
Як ми переконались, описаний вище метод Стеффенсена в класі сильно випуклих функцій має квадратичну швидкість збіжності Розглянутий метод Стеффенсена, можна використовувати для безумовної мінімізації функцій багатьох змінних. Але недоліком методу (2) являється велика необхідність в досить великій близкості початкової точки u0.