Смешанная задача для уравнения гиперболического типа

Здесь u_ij - приближенное значение функции u(x,t) в узле (x_i,t_j).

Полагая, что l = t / h , получаем трехслойную разностную схему

u_i,j+1 = 2(1- l ² )u_i,j + l ² (u_i+1,j- u_i-1,j) - u_i,j-1 , i = 1,2 ... n. (5)

Для простоты в данной лабораторной работе заданы нулевые граничные условия, т.е. m ₁(t) º 0, m ₂(t) º 0. Значит, в схеме (5) u_0,j= 0, u_nj=0 для всех j. Схема (5) называется трехслойной на трех временных слоях с номерами j-1, j , j+1. Схема (5) явная, т.е. позволяет в явном виде выразить u_i,j через значения u с предыдущих двух слоев.

Численное решение задачи состоит в вычислении приближенных значений u_i,j решения u(x,t) в узлах (x_i,t_j) при i =1, ... n, j=1,2, ... ,m . Алгоритм решения основан на том, что решение на каждом следующем слое ( j = 2,3,4, ... n) можно получить пересчетом решений с двух предыдущих слоев ( j=0,1,2, ... , n-1) по формуле (5). На нулевом временном слое (j=0) решение известно из начального условия u_i0 = f(x_i).

Для вычисления решения на первом слое (j=1) в данной лабораторной работе принят простейший способ, состоящий в том, что если положить ¶ u(x,0)/ ¶ t ^» ( u( x, t ) - u(x,0) )/ t (6) , то u_i1=u_i0+ + t (x_i), i=1,2, ... n. Теперь для вычисления решений на следующих слоях можно применять формулу (5). Решение на каждом следующем слое получается пересчетом решений с двух предыдущих слоев по формуле (5).

Описанная выше схема аппроксимирует задачу с точностью до О( t +h²). Невысокий порядок аппроксимации по t объясняется использованием слишком грубой аппроксимации для производной по е в формуле (6).

Схема устойчива, если выполнено условие Куранта t < h. Это означает, что малые погрешности, возникающие, например, при вычислении решения на первом слое, не будут неограниченно возрастать при переходе к каждому новому временному слою. При выполнении условий Куранта схема обладает равномерной сходимостью, т.е. при h ® 0 решение разностной задачи равномерно стремится к регшению исходной смешанной задачи.

Недостаток схемы в том, что как только выбраная величина шага сетки h в направлении x , появляется ограничение на величину шага t по переменной t . Если необходимо произвести вычисление для большого значения величины T , то может потребоваться большое количество шагов по переменной t. Указанный гнедостаток характерен для всех явных разностных схем.

Для оценки погрешности решения обычно прибегают к методам сгущения сетки.

Для решения смешанной задачи для волнового уравнения по явной разностной схеме (5) предназначена часть программы, обозначенная Subroutine GIP3 Begn ... End . Данная подпрограмма вычисляет решение на каждом слое по значениям решения с двух предыдущих слоев.

Входные параметры :

hx - шаг сетки h по переменной х;

ht - шаг сетки t по переменной t;

k - количество узлов сетки по x, a = hn;

u1 - массив из k действительных чисел, содержащий значение решений на ( j - 1 ) временном слое, j = 1, 2, ... ;

u2 - массив из n действительных чисел, содержащий значение решений на j - м временном слое, j = 1, 2, ... ;

u3 - рабочий массив из k действительных чисел.

Выходные параметры :

u1 - массив из n действительных чисел, содержащий значение решения из j - м временном слое, j = 1, 2, ... ;

u2 - массив из n действительных чисел, содержащий значение решения из ( j +1) - м временном слое, j = 1, 2, ... .

К части программы, обозначенной как Subroutine GIP3 Begin ... End происходит циклическое обращение, пеоред первым обращением к программе элементам массива u2 присваиваются начальные значения, а элементам массива u1 - значения на решения на первом слое, вычислинные по формулам (6). При выходе из подпрограммы GIP3 в массиве u2 находится значение решения на новом временном слое, а в массиве u1 - значение решения на предыдущем слое.

Порядок работы программы:

1) описание массивов u1, u2, u3;

2) присвоение фактических значений параметрам n, hx, ht, облюдая условие Куранта;

3) присвоение начального значения решения элементам массива и вычисленное по формулам (6) значение решения на первом слое;

4) обращение к GIP3 в цикле k-1 раз, если требуется найти решение на k-м слое ( k ³ 2 ).

Пример:

Решить задачу о колебании струны единичной длины с закрепленными концами, начальное положение которой изображено на рисунке. Начальные скорости равны нулю. Вычисления выполнить с шагом h по x, равным 0.1, с шагом t по t, равным 0.05, провести вычисления для 16 временных слоев с печатью результатов на каждом слое. Таким образом, задача имеет вид

( ¶ ² u/ ¶ t²) = ( ¶ ² u/ ¶ x ²) , x Î [ 0 , 1 ] , t Î [ 0 , T ] ,

u ( x , 0 ) = f (x) , x Î [ 0 , a ], ¶ u(x,0)/ ¶ t = g(x) , x Î [ 0 , a ],

u ( 0 , t ) = 0, u ( 1 , t ) = 0, t Î [ 0 , 0.8 ],

æ 2x , x Î [ 0 , 0.5 ] ,

f(x) = í g( x ) = 0

î 2 - 2x , x Î [ 0.5 , 1 ] ,

{ ************************************************************* } 
{ Приложение 3  ( выполнения лабораторной работы. Вариант 12)   } 
{ ------------                                                  } 
{  Программа решения смешанной задачи для уравнения гиперболи-  } 
{  ческого типа методом сеток.                                  } 
{  Выполнил студент гр. МС-2136    Осинцев А.В.                 } 
{ ************************************************************* } 
 
Program Laboratornaya_rabota_43_variant_12; 
 
Const 
 
  hx = 0.1 ;     { Шаг по x - hx } 
  ht = 0.05 ;    { Шаг по t - ht } 
  n  = 11 ;      { Количество узлов } 
 
Function f(x : Real) : Real; { Данная функция              } 
                             { вычисляющая решение при t=0 } 
Begin 
 
  f := sin(pi * x) * cos(x); 
 
End; 
 
Function g(x : Real) : Real; { Данная функция                          } 
                             { вычисляющая производную решения при t=0 } 
Begin 
  g := 0; 
End; 
 
Var 
  xp            : Array[1..n] of Real; 
  i,j,n1        : Word; 
  x,t,a1,b1     : Real; 
  u1,u2,u3      : Array[1..n] of Real; 
 
Begin 
  n1 := n; 
  WriteLn('Приложение 4'); 
  WriteLn('------------'); 
  WriteLn('Результат, полученный при вычислении программы :'); 
  WriteLn; 
  xp[1] := 0; 
  xp[n] := 1; 
 
  For i := 2 to ( n - 1 ) do 
    Begin 
      x := (i-1) * hx; 
      xp[i] := x; 
      u1[i] := f(x);                 { u(x,0) на 0 слое } 
      u2[i] := u1[i] + ht * g(x);    { u(x,ht) на 1 слое } 
    End; 
 
  {   ///   Задание граничных условий  \\\   } 
 
  u1[1] := 0 ;   {  u(0,0)   } 
  u1[n] := 0 ;   {  u(1,0)   } 
  u2[1] := 0 ;   {  u(0,ht)  } 
  u2[n] := 0 ;   {  u(1,ht)  } 
  u3[1] := 0 ;   {  u(0,2ht) } 
  u3[n] := 0 ;   {  u(1,2ht) } 
 
  {   ///   Печать заголовка   \\\   } 
  Write('       '); 
  For i := 1 to n do Write(' x=', xp[i]:1:1); 
  WriteLn; 
 
  t := 0; 
  {   ///   Печать решения на нулевом слое   \\\   } 
  Write('t=',t:2:2,' '); 
  For i := 1 to n do 
    If u1[i] >= 0 then Write(' ',u1[i]:3:3) else Write(u1[i]:3:3) ; 
 
  t := t + ht; 
  {   ///   Печать решения на первом слое   \\\   } 
  WriteLn; 
  Write('t=',t:2:2,' '); 
  For i := 1 to n do 
    If u2[i] >= 0 then Write(' ',u2[i]:3:3) else Write(u2[i]:3:3); 
 
  For j := 1 to 15 do 
    Begin 
      {Subroutine GIP3 Begin} 
      n1 := n1-1; 
      {Вычисление параметра сетки для проверки условия Куранта} 
      a1 := ht/hx; 
      if a1 > 1 then WriteLn('Нарушено условие Куранта') else 
        Begin 
          b1 := a1 * a1; 
          a1 := 2 * ( 1 - b1); 
          {Вычисление решения на очередном слое} 
          For i := 2 to n do u3[i] := a1*u2[i] + b1 * (u2[i+1] + 
                             u2[i-1]) - u1[i]; 
          For i := 2 to n do 
            Begin 
              u1[i] := u2[i]; 
              u2[i] := u3[i] 
            End; 
 
        End; 
      u1[n] := 0; 
      u2[n] := 0; 
      u3[n] := 0; 
      {Subroutine GIP3 End} 
    t := t + ht; 
    WriteLn; 
    Write('t=',t:2:2,' '); 
    For i := 1 to n do 
      {Вывод результатов} 
      If u2[i] >= 0 then Write(' ',u2[i]:3:3) else Write(u2[i]:3:3); 
 
    End; 
 
  WriteLn; 
  WriteLn; 
End.     

Дата добавления: 30.10.2000