197-577-902

Лекции по информатике - Метод итераций (повторений)

Пусть М – некоторое множество (значения должны быть допустимы ), Р подмножество М, где Р – множество решений. При этом . Пусть Т некоторое преобразование из М\Р в М.

Метод итераций (дословно – повторений) состоит в том, что строится некоторое преобразование (отображение) и это преобразование последовательно применяется, начиная с какого-то :

до тех пор, пока мы впервые не получим x_i такое, что (естественно, должна быть уверенность, что рано или поздно такое x_i найдется).

Напишем программу, находящую х. Пусть (некоторая точка М) – начальное значение.

x:= х0;
while not p(x) do
begin
  x:= T(x);
end;

Для корректности данной программы необходимо и достаточно доказать, что цикл завершится за конечное число шагов.

Иллюстрация работы программы:

Иллюстрация работы программы

При этом множества М и Р могут быть выбраны различными способами.

Вычисление интеграла Римана

Для вычисления интеграла необходимо разбить отрезок [a,b] на р частей. При этом , где I_p– квадратурная форма. Если увеличить количество разбиений в к раз, то получим . При этом если подынтегральная функция достаточно гладкая, то , где s – степень точности формулы. Данная формула оценки погрешности через старший член лежит в основе метода Рунге.

При этом величина будет точнее, чем I_kp. Она носит название первого шага метода Ромберга.

Рассмотрим программу вычисления интеграла.

{$N+}
type
  _Real = single;
  fr2r = function (x: _Real):_Real;


function Integral (a, b: _Real; f: fr2r; eps: Real):Real;
{ a, b – пределы интегрирования, f – достаточно гладкая
  подынтегральная функция; 
  eps>0 – критерий Рунге.
  Функция возвращает приближенное значение интеграла Римана }
var
  Iprev, Ipost, R: _Real;
  p: word;
  h: _Real;
const
  p0 = 16;
  k = 2; {или 3}
  gamma = ... { k_s+1 - 1}

{ далее могут следовать другие необходимые константы }
begin
  if a = b then
  begin
    Integral:= 0.0;
    exit;
  end;
  p:= p0;
  h:= (b-a) / p;

{ Вычисление Ipostдля числа разбиений p, а также,
  возможно, его частей }
  R:=1e30;
  while abs(R) > eps do
  begin
    Iprev:= Ipost;
    p:= p*k;
    h:= h / k; { p-иногда не обязателен }
  { Вычисляем Ipost для pразбиений, и, возможно, его части }
    R:= (Ipost - Iprev) / gamma;
  end;
  Integral:= Ipost + R;
end;

Замечания:

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

Вычисление функции Бесселя .

С другой стороны при достаточно больших N.

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

function J0(x, eps : _Real):_Real;
{ x – аргумент функции, eps>0 – погрешность округления.
  Функция возвращает приближенное значение функции Бесселя }
var
  k: word;
  a, s: _Real;
begin
  k:= 0;
  a:=1.0;
  s:= 1.0;
  while abs(a) > eps do
  begin
    k:= k + 1;
    a:= -a * sqr (x / 2.0) / sqr (k);
    s:= s + a;
  end;
  j0 := s;
end;

При этом k лучше было бы хранить как вещественное число.

Для увеличения скорости работы программы удобно затабулировать квадрат знаменателя. Запишем разность знаменателей на следующих один за другим шагах цикла. Z_k₊₁ = (k+1)²-k² = 2k + 1, в свою очередь, табулируя полученную линейную функцию, получаем U_k₊₁= Z_k₊₁ –Z_k = 2. Таким образом нахождение k² сводится к сложению представленному следующей таблицей:

k	0	1	2	3	4
k²	0	1	4	9	16
z	-1	1	3	5	7
u	2	2	2	2	2

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

type
  _Real = single;

function J0 (x, eps : _Real) : _Real;
{ x – значение аргумента, eps>0 – погрешность округления.
  Функция возвращает приближенное значение функции Бесселя }
const
  u : _Real = 2.0;
var
  s, a, z, z0, r : _Real;
begin
  a:= 1.0;
  s:= a;
  z:= -1.0;
  z0:= 0.0;
  r:= -sqr (x/2.0);
  while abs (a) > eps do
  begin
    z:= z + u;
    z0 := z0 + z;
    a:= a * r / z0;
    s:= s+a;
  end;
  J0:= s;
end;

Математические свойства данной программы будут те же, т.к. используется тот же математический аппарат. Функция Бесселя, вычисляемая в данной программе характеризуется следующими графиками :

, где t – длина мантиссы, b - основание системы счисления принятой в компьютере.

Таким образом, не следует брать слишком маленькие значения eps. Минимальное значение eps зависит от х и возрастает с ростом модуля х.

Из всего вышесказанного можно сделать вывод, что метод рядов Тейлора хотя и является хорошим математическим аппаратом для исследований, но не приспособлен для вычислений, т.к. область его применения ограничивается небольшой окрестностью точки разложения (в данном случае точки 0).