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

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

Рассмотрим смешанную задачу для волнового уравнения ( 2 u/ t2) = c 2 * ( 2u/ x2) (1). Задача состоит в отыскании функции u(x,t) удовлетворяющей данному уравнению при 0 < x < a, 0 < t T, начальным условиям u(x,0) = f(x), u(x,0)/ t = g(x) , 0 x a и нулевыми краевыми условиями u(0,t) = u(1,t)=0.

Так как замена переменных t  ct приводит уравнение (1) к виду ( 2 u/ t2) = ( 2u/ x2), то в дальнейшем будем считать с = 1.

Для построения разностной схемы решения задачи строим в области D = {(x,t) | 0 x a, 0 t T } сетку x>i> = ih, i=0,1 ... n , a = h * n, t>j> = j*  , j = 0,1 ... , m,  m = T и аппроксимируем уравнение (1) в каждом внутреннем узле сетки на шаблоне типа “крест”.

t

T

j+1

j

j-1

0 i-1 i i+1


Используя для аппроксимации частных производных центральные разностные производные, получаем следующую разностную аппроксимацию уравнения (1) .

u>i,j+1> - 2u>ij> + u>i,j-1> u>i+1,,j> - 2u>ij> + u>i-1, j>


2 h2


(4)

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

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

u>i,j+1> = 2(1-  2 )u>i,j> +  2 (u>i+1,j>- u>i-1,j>) - u>i,j-1> , i = 1,2 ... n. (5)

Для простоты в данной лабораторной работе заданы нулевые граничные условия, т.е.  >1>(t)  0,  >2>(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,  ) - u(x,0) )/  (6) , то u>i1>=u>i0>+ +  (x>i>), i=1,2, ... n. Теперь для вычисления решений на следующих слоях можно применять формулу (5). Решение на каждом следующем слое получается пересчетом решений с двух предыдущих слоев по формуле (5).

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

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

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

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

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

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

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

ht - шаг сетки  по переменной 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, ... .

К части программы, обозначенной как sub>routine 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 ).

Пример:

1

0.5 0.5

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

(  2 u/  t2) = (  2 u/  x 2) , 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] ,

Строим сетку из 11 узлов по x и выполняем вычисления для 16 слоев по t. Программа, и результаты вычисления приведены далее.

Приложение 1

(пример выполнения лабораторной работы)

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

Program Laboratornaya_rabota_43;

Const

hx = 0.1 ; { Шаг по x - hx }

ht = 0.05 ; { Шаг по t - ht }

n = 11 ; { Количество узлов }

Function f(x : Real) : Real; { Данная функция }

{ вычисляющая решение при t=0 }

Begin

If x <= 0.5 then

f := 2 * x

else

f := 2 - 2 * 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('Приложение 2');

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

{sub>routine 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;

{sub>routine 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.

Приложение 3

( выполнения лабораторной работы. Вариант 11)

Program Laboratornaya_rabota_43_variant_11;

Const

hx = 0.1 ; { Шаг по x - hx }

ht = 0.05 ; { Шаг по t - ht }

n = 11 ; { Количество узлов }

Function f(x : Real) : Real; { Данная функция }

{ вычисляющая решение при t=0 }

Begin

f := x * ( x * x - 1 );

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

{sub>routine 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;

{sub>routine 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.

(выполнения лабораторной работы. Вариант 20)

Program Laboratornaya_rabota_43_variant_20;

Const

hx = 0.1 ; { Шаг по x - hx }

ht = 0.05 ; { Шаг по t - ht }

n = 11 ; { Количество узлов }

Function f(x : Real) : Real; { Данная функция }

{ вычисляющая решение при t=0 }

Begin

f := 10 * x * ( x * x * x - 1 );

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

{sub>routine 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;

{sub>routine 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.

( выполнения лабораторной работы. Вариант 14)

Program Laboratornaya_rabota_43_variant_14;

Const

hx = 0.1 ; { Шаг по x - hx }

ht = 0.05 ; { Шаг по t - ht }

n = 11 ; { Количество узлов }

Function f(x : Real) : Real; { Данная функция }

{ вычисляющая решение при t=0 }

Begin

f := x * sin( 2 * (x - 1) );

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

{sub>routine 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;

{sub>routine 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.

( выполнения лабораторной работы. Вариант 13)

Program Laboratornaya_rabota_43_variant_13;

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 * 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

{sub>routine 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;

{sub>routine 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.

( выполнения лабораторной работы. Вариант 12)

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

{sub>routine 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;

{sub>routine 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.