Полином Гира (экстраполяция методом Гира)

{****************************************************************************}

{                                                                            }

{ Вариант №1.23   Экстраполяция методом Гира         гр. ##-###       Ф.И.О. }

{                                                                            }

{           P^2                                     P = 50-200 [Ом],         }

{   f(w) = ----- * ( 1 + Q^2 * (2*dw/wp))^(-1/2),   Rn= 0.1-5  [Ом],         }

{           Rn                                      Q = 50-1000,             }

{                                                   wp- частота резонанса,   }

{                                                   dw- шаг,                 }

{                                                                            }

{****************************************************************************}

{                              Аппроксимация

      Аппроксимацией   называется  замена   одних  математических  объектов

   (чисел, функций)  другими,  более  простыми  в том  или  ином  смысле  и

   близкими к исходным.

      При  решении  задачи   аппроксимации  предполагается,   что  известны

   значнения  аппроксимируемой  функции  f(x)  (снятые эксперементально или

   расчитанные по формуле)  в  узловых  точках  (обычно их называют  узлами

   интерполяции) X1, X2, ..., Xn.  Необходимо найти функцию P(x), с помощью

   которой можно  было бы определить  значения исходной функции  для других

   (отличных от узлов)  значений аргумента  x  с некоторой достаточно малой

   погрешностью.

      Различают две задачи аппроксимации: интерполяцию и экстраполяцию.

      Интерполяция  -  это отыскание промежуточных  значений функции внутри

   заданного интервала [X1, Xк]  для значений  x  не совпадающих с узловыми

   точками.

      Экстраполяция - это распространение выводов, полученных из наблюдения

   над одной  частью явления  на другую ее часть,  то есть  это определение

   значений функции вне заданного интервала [Xн, Xк]. }

PROGRAM GIR;

{$N-}

USES Crt;

CONST

   Max_kol = 25;                       {Максимальное количество узлов}

   Tochka  = 10;                       {Количество знаков после запятой}

TYPE

   _Real = Real;

   _INTEGER = Byte;

   _Mas = array [1..Max_kol] of _Real;

VAR

   H, F_up, F_down, eps_pol, eps_tek   : _Real;

   Xpol, Ypol, V                       : _Mas;

   n, m                                : _INTEGER;     {n-количество узлов}

   Pr                                  : Char;         {переключатель}

{----------------------------------------------------------------------------}

{Расчет табличных данных для заданной функции}

PROCEDURE TABL_FUNC(Var Xpol, Ypol: _Mas; Var n: _INTEGER; Var dw, eps_pol: _Real);

VAR

   P, Rn, Q, wp, wn   : _Real;

   I                  : _INTEGER;            {счетчик}

BEGIN

   Write('Сопротивление P ............. = '); ReadLn(P);

   Write('Сопротивление нагрузки Rn ... = '); ReadLn(Rn);

   Write('Параметр Q .................. = '); ReadLn(Q);

   Write('Начальна частота wp ......... = '); ReadLn(wn);

   Write('Шаг изменения частота dw .... = '); ReadLn(dw);

   Write('Количество узлов (< ', Max_kol:2, ') ..... = '); ReadLn(n);

   Write('Точность eps_pol ............ = '); ReadLn(eps_pol);

   WriteLn;

  {расчет функции в узлах}

   wp:= wn;

   FOR I:= 1 TO (n+1) DO

    begin

       Xpol[I]:= wp;

       Ypol[I]:= Sqr(P) / (Rn * Sqrt(1 + Sqr(Q) * Sqr(2*dw/wp)));

       wp:= wp + dw;

    end;

END;

{--------------------------< Конец ПП TABL_FUNC >----------------------------}

{Ввод значений аппроксимируемой функции}

PROCEDURE VVOD_TABL(Var Xpol, Ypol: _Mas; Var n: _INTEGER; Var dw, eps_pol: _Real);

VAR

   wp, wn   : _Real;

   I        : _INTEGER;                      {счетчик}

BEGIN

   Write('Начальный узел ............ = '); ReadLn(wn);

   Write('Шаг изменения узла ........ = '); ReadLn(dw);

   Write('Количество узлов (< ', Max_kol:2, ') ... = '); ReadLn(n);

   Write('Погрешность eps .............. = '); ReadLn(eps_pol);

   WriteLn;

  {ввод значений функции в узлах}

   WriteLn('  Введите значение функции в узлах:');

   wp:= wn;

   FOR I:= 1 TO n DO

    begin

       Xpol[I]:= wp;

       Write('Y(', wp:Tochka, ')= '); ReadLn(Ypol[I]);

       wp:= wp + dw;

    end;

  {задание (n+1)-узла}

   Xpol[n+1]:= wp;

END;

{---------------------------< Конец ПП VVOD_TABL >---------------------------}

{Расчет коэффициентов полинома}

PROCEDURE Kff_pol(X: _Mas; H: _Real; n, m: _INTEGER; Var V: _Mas);

TYPE

   _Matric = array [1..Max_kol] of _Mas;

   _Vector = array [1..Max_kol] of Byte;

VAR

   G      : _Matric;

   P      : _Vector;

   I, J   : _INTEGER;                        {счетчики}

{----------------------------------------------------------------------------}

{рашение матричной системы уравнений, методом "LU-разложений"}

PROCEDURE GVP(G: _Matric; P: _Vector; m: _INTEGER; Var V: _Mas);

VAR

   K, I, J   : _INTEGER;                         {счетчики}

   L, U      : _Matric;                          {вспомогательные матрицы}

   SOS       : _Mas;                             {вспомогательный вектор}

   alf, S    : _Real;                            {вспомогательные переменные}

BEGIN

  {обнуление элементов вспомогательных матриц L и U}

   FOR I:= 1 TO m DO

   For J:= 1 To m Do

    begin

       L[I, J]:= 0;                              { L[I, J] = 0, где I=1..m, J=1..m   }

       U[I, J]:= 0;                              { U[I, J] = 0, где I=1..m, J=1..m   }

    end;

    {разложение матрицы G на матрицу L и U}

     FOR K:= 1 TO m DO

      begin

        {формирование матрицы L}

         For I:= K To m Do

          begin

            {накопление суммы элементов в S}

             S:= 0;

             for J:= 1 to (K - 1) do

              begin                              {    k-1                                      }

                 alf:= L[I, J] * U[J, K];        { S = E (L[I, J] * U[J, K]),                  }

                 S:= S + alf;                    {    j=1                                      }

              end;

             L[I, K]:= G[I, K] - S;              { L[I, K] = G[I, K] - S, где k=1..m, i=k..m   }

          end;

        {формирование матрицы U}

         U[K, K]:= 1;

         For I:= (K + 1) To m Do

          begin

            {накопление суммы элементов в S}

             S:= 0;

             for J:= 1 to (K - 1) do

              begin                              {    k-1                                             }

                 alf:= L[K, J] * U[J, I];        { S = E (L[K, J] * U[J, I]),                         }

                 S:= S + alf;                    {    j=1                                             }

              end;

                                                 {            G[I, K] - S                                           }

             alf:= G[I, K] - S;                  { U[K, I] = -------------,  где k=1..m, i=(k+1)..m   }

             U[K, I]:= alf/L[K, K];              {              L[K, K]                                             }

          end;

      end;

    {вычисление вспомогательного вектора SOS}

     FOR K:= 1 TO m DO

      begin

        {накопление суммы элементов в S}

         S:= 0;

         For I:= 1 To (K - 1) Do

          begin                                  {    k-1                            }

             alf:= L[K, I] * SOS[I];             { S = E (L[K, I] * SOS[I]),         }

             S:= S + alf;                        {    i=1                            }

          end;

                                                 {           P[K] - S                }

         alf:= P[K] - S;                         { SOS[K] = ----------, где k=1..m   }

         SOS[K]:= alf/L[K, K];                   {            L[K, K]                }

      end;

    {вычисление вектора V}

     FOR K:= m DOWNTO 1 DO

      begin

        {накопление суммы элементов в S}

         S:= 0;

         For I:= (K + 1) To m Do

          begin                                  {       m                         }

             alf:= U[K, I] * V[I];               { S =   E  (U[K, I] * V[I]),      }

             S:= S + alf;                        {     i=k+1                       }

          end;

         V[K]:= SOS[K] - S;                      { V[K] = SOS[K] - S, где k=1..m   }

      end;

  END;

{------------------------------< конец ПП GVP >------------------------------}

BEGIN

 {формирование матрицы G}

  FOR I:= 1 TO m DO                                        {           | X[n+1] - X[n - (m-1)] | (m-1)   }

   For J:= 1 To m Do                                       { G[I, J] = |-----------------------|         }

   G[I, J]:= exp( (I-1)*Ln((X[n+1] - X[n-(J-1)]) / H) );   {           |            H          |         }

 {формирование вектора P = [1, 0, 0, 0, ..., 0]}

  P[1]:= 1;

  FOR I:= 2 TO m DO P[I]:= 0;

 {обращение к ПП GVP}

  GVP(G, P, m, V);

END;

{---------------------------< Конец ПП Kff_pol >-----------------------------}

{Расчет полинома}

FUNCTION F_pol(V, Y: _Mas; n, m: _INTEGER): _Real;

VAR

   I: _INTEGER;                              {счетчик}

   S: _Real;                                 {сумматор}

BEGIN

   S:= 0;

   FOR I:= 1 TO m DO S:= S + V[I]*Y[n+1-I];

   F_pol:= S;

END;

{-----------------------------< Конец ПП F_pol >-----------------------------}

{Вывод значения функции в (n+1) узле}

PROCEDURE VIVOD1(Xpol, Ypol, eps_tek: _Real; n, m: _INTEGER);

BEGIN

   WriteLn;

   WriteLn;

   WriteLn('Значение функции в (n+1) узле:  Y(', Xpol:Tochka, ')=', Ypol:Tochka);

   WriteLn('Количество используемых узлов', m:3, ' из', n:3);

   WriteLn('Погрешность расчета =', eps_tek:Tochka);

END;

{-----------------------------< Конец ПП VIVOD1 >----------------------------}

{Вывод значения функции в (n+1) узле для тестовой функции}

PROCEDURE VIVOD2(Xpol, Ypol, eps_tek, Y: _Real; n, m: _INTEGER);

VAR

   eps: _Real;

BEGIN                                            {               Y_полином - Y_точное  }

   eps:= Abs((Ypol - Y)/Y);                      { EPS_точное = ---------------------- }

                                                 {                     Y_точное        }

   WriteLn;

   WriteLn('  Значение функции в узле X(n+1)=', Xpol:Tochka , ' :');

   WriteLn('Y полинома ... =', Ypol:Tochka);

   WriteLn('Y точное ..... =', Y:Tochka);

   WriteLn('  Погрешность расчета :');

   WriteLn('текущая ...... =', eps_tek:Tochka);

   WriteLn('точная ....... =', eps:Tochka);

   WriteLn;

   WriteLn('Количество используемых узлов', m:3, ' из', n:3);

END;

{-----------------------------< Конец ПП VIVOD2 >----------------------------}

{основной блок программы}

BEGIN

   ClrScr;

   WriteLn('  ExtraPol_GIRA 1.00          Copyright (c) 1999-00                      DGur.');

   WriteLn;

   WriteLn;

  {Ввод исходных данных}

   Write('Тестовый расчет? [Y/N]'); ReadLn(Pr);

   WriteLn;

  {Ввод значений функции в узлах}

   IF (UpCase(Pr) = 'Y')

     {обращение к ПП TABL_FUNC}

      THEN TABL_FUNC(Xpol, Ypol, n, H, eps_pol)

     {обращение к ПП VVOD_TABL}

      ELSE VVOD_TABL(Xpol, Ypol, n, H, eps_pol);

  {итерационный процесс}

   m:= 2;

   F_down:= 0;

   REPEAT

     {обращение к ПП Kff_pol}

      Kff_pol(Xpol, H, n, m, V);

     {обращение к ПФ F_pol}

      F_up:= F_pol(V, Ypol, n, m);

     {точность расчета}

      IF Not(F_up = 0) THEN eps_tek:= Abs((F_up - F_down)/F_up);

      F_down:= F_up;

      Inc(m);

   UNTIL (eps_tek < eps_pol) or (m > n);

  {вывод полученных результатов}

   IF (UpCase(Pr) = 'Y')

     {обращение к ПП VIVOD2}

      THEN VIVOD2(Xpol[n+1], F_up, eps_tek, Ypol[n+1], n, (m-1))

     {обращение к ПП VIVOD1}

      ELSE VIVOD1(Xpol[n+1], F_up, eps_tek, n, (m-1));

   ReadKey;

END.                                             {Конец основного блока}