Вычислительные методы: разностная схема и метод прогонки (Delphi)

program Lab7_3;

{$APPTYPE CONSOLE}

{ Описание задания по вычислительным методам.

 Дана разностная схема
   A[I]*Y[I+1]+B[I]*Y[I]+C[I]*Y[I-1] = 12*H^2*F(x[I]) + H^4*F''(x[I]), где
      I = 1,...n-1
      F'' - вторая производная
   Y[0] = 0
   Y[n] = 0
   H = 1/n
   x[I] = I * H, где I = 1,...n-1

   A[I] = H^2 - 12, I = 1,...n-1
   B[I] = -(10*H^2 + 24), I = 1,...n-1
   C[I] = A[I], I = 1,...n-1

 аппроксимирующая задачу -d^2 u/dx^2 + u = F(x), где
     d - знак дифференциала
     x принадлежит отрезку [0, 1]
     u(0) = 0
     u(1) = 0
 с четвертым порядком.


 Полагаем u(x) = Sin(2*Pi*x).

 Необходимо решить разностную схему методом прогонки и сравнить
 полученные значения с истинными
}

uses
  SysUtils;

type
  TFloat = Extended;
  TFloatArray = array of TFloat;

function Progonka(A, B, C, F: TFloatArray; var Y: TFloatArray): Boolean;
// Решение системы уравнений с трехдиагональной матрицей
// методом правой прогонки.
//
// Обозначим за X[i] - Х с индексом i.
//   A[i] * Y[I+1] - B[I] * Y[I] + C[I] * Y[I-1] = F[I], где I от 1 до N-1,
//                                                       A[I] <> 0
//                                                       C[I] <> 0
//   Y[0] = A[0] * Y[1] + F[0]
//   Y[N] = C[N] * Y[N-1] + F[N]
//
// Т.о. данная система уравнений на Y-итые образует трехдиагональную матрицу.
// Решить данную систему можно используя следующие формулы:
//   K[0]:= A[0]
//   L[0]:= F[0]
//
//   K[I]:= A[I] / (B[I] - K[I-1]*C[I]), где I от 1 до N-1
//   L[I]:= (C[I]*L[I-1] - F[I]) / (B[I} - K[I-1]*C[I]), где I от 1 до N-1
//
//   Y[N]:= (F[N] + C[N]*L[N-1]) / (1 - C[N]*K[N-1])
//   Y[I]:= K[I]*Y[I+1] + L[I], где I от N-1 до 0
var
  n, I: Integer;
  K, L: TFloatArray;
begin
  if Length(A) < 3 then
  begin
    Result:= False;
    Exit;
  end;

  n:= Length(A) - 1;
  if (Length(B) <> n + 1) or (Length(C) <> n + 1) or (Length(F) <> n + 1) then
  begin
    Result:= False;
    Exit;
  end;

  try
    SetLength(K, n);
    SetLength(L, n);
    SetLength(Y, n + 1);

    K[0]:= A[0];
    L[0]:= F[0];

    for I:= 1 to n - 1 do
    begin
      K[I]:= A[I] / (B[I] - K[I-1] * C[I]);
      L[I]:= (C[I] * L[I-1] - F[I]) / (B[I] - K[I-1] * C[I]);
    end;

    Y[n]:= (F[n] + C[n] * L[n-1]) / (1 - C[n] * K[n-1]);
    for I:= n - 1 downto 0 do
      Y[I]:= K[I] * Y[I+1] + L[I];

    Result:= True;
  except
    Result:= False;
  end;
end;

procedure DoWork(A, B: TFloat; n: Integer);
// Функция выполняет задание
//   A, B - концы отрезка
//   n - кол-во интервалов на отрезке
var
  _a, _b, _c: TFloatArray; // коэффициенты в левой части
  _f: TFloatArray;         // значения правой части в узлах
  Y: TFloatArray;          // приближенное решение полученное
                           // методом прогонки
  H: TFloat;               // шаг
  I: Integer;
  Err: TFloat;
begin
  if A >= B then
    Exit;

  if n < 3 then
    Exit;

  H:= Abs(B-A) / n;

  try
    SetLength(_a, n + 1);
    SetLength(_b, n + 1);
    SetLength(_c, n + 1);
    SetLength(_f, n + 1);

    for I:= 0 to n do
    begin

      // Задаем коэффициенты разностной схемы

      if I = 0 then
        _a[I]:= 0 // следует из нашей разностной схемы
      else
        _a[I]:= H*H-12;

      _b[I]:= -(10*H*H+24);

      if I = n then
        _c[I]:= 0 // следует из нашей разностной схемы
      else
        _c[I]:= H*H-12;


      // Задаем правую часть в разностной схеме

      if I = 0 then
        _f[I]:= 0 // следует из нашей разностной схемы
      else if I = n then
        _f[I]:= 0 // следует из нашей разностной схемы
      else
        _f[I]:= Sqr(H)*(1+4*Sqr(Pi))*Sin(2*Pi*(A+I*H))*(12 - Sqr(H*2*Pi));
    end;

    // Решаем полученную систему методом прогонки

    if Progonka(_a, _b, _c, _f, Y) then
    begin
      Err:= 0;
      for I:= 0 to n do
      begin

        // Сравнимаем полученное решение с истинным решением

        Write('Y[', I, '] = ', Y[I]);
        Write('':10, 'U[', I, '] = ', Sin(2*Pi*(A+I*H)));

        // Вычисляем точность, как максимум модуля разности
        // между полученным (приближенным) решением и истинным решением

        if Abs(Y[I] - Sin(2*Pi*(A+I*H))) > Err then
          Err:= Abs(Y[I] - Sin(2*Pi*(A+I*H)));

        WriteLn;
      end;

      WriteLn('Mistake: ', Err);

    end else
     WriteLn(ErrOutput, 'Unexpected error!!!');
  except
    on E: Exception do
      WriteLn(ErrOutput, E.Message);
  end;
end;

begin
  DoWork(0, 1, 10);
  ReadLn;
end.
 

Слава Антонов © 2002 — August 13, 2008

197-577-902