Смешанная задача для уравнения гиперболического типа
МИНИСТЕРСТВО ОБЩЕГО И ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ Р.Ф.
КУРГАНСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
Кафедра прикладной и высшей математики
Лабораторная работа № 43
на тему:
Решение смешанной задачи
для уравнения
гиперболического типа
методом сеток
Группа М-2136
Выполнил студент _______________________
Проверил преподаватель Воронова Лилия Ивановна
Курган 1998
Рассмотрим смешанную задачу
для волнового уравнения ( ¶ 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 } сетку xi = ih, i=0,1 ... n , a = h * n, tj
= j* ttt , j = 0,1 ... , m, t m = T и
аппроксимируем уравнение (1) в каждом внутреннем узле сетки на шаблоне типа
“крест”.
Используя для аппроксимации
частных производных центральные разностные производные, получаем следующую
разностную аппроксимацию уравнения (1) .
ui,j+1 - 2uij + ui,j-1
ui+1,,j - 2uij + ui-1, j
t 2 h2
|
(4)
Здесь uij -
приближенное значение функции u(x,t) в узле (xi,tj).
Полагая, что l = t / h , получаем трехслойную разностную схему
ui,j+1 = 2(1- l 2 )ui,j + l 2 (ui+1,j-
ui-1,j) - ui,j-1 , i = 1,2 ... n. (5)
Для простоты в данной
лабораторной работе заданы нулевые граничные условия, т.е. m 1(t) º 0, m 2(t) º 0. Значит, в схеме (5) u0,j=
0, unj=0 для всех j. Схема (5) называется трехслойной на трех
временных слоях с номерами j-1, j , j+1. Схема (5) явная, т.е. позволяет в
явном виде выразить ui,j через значения u с предыдущих двух слоев.
Численное решение задачи
состоит в вычислении приближенных значений ui,j решения u(x,t) в
узлах (xi,tj) при i =1, ... n, j=1,2, ... ,m . Алгоритм
решения основан на том, что решение на каждом следующем слое ( j = 2,3,4, ...
n) можно получить пересчетом решений с двух предыдущих слоев ( j=0,1,2, ... ,
n-1) по формуле (5). На нулевом временном слое (j=0) решение известно из
начального условия ui0 = f(xi).
Для вычисления решения на
первом слое (j=1) в данной лабораторной работе принят простейший способ,
состоящий в том, что если положить ¶ u(x,0)/ ¶ t » ( u( x, t ) - u(x,0) )/ t (6) , то ui1=ui0+
+ t (xi), i=1,2, ... n. Теперь для вычисления
решений на следующих слоях можно применять формулу (5). Решение на каждом
следующем слое получается пересчетом решений с двух предыдущих слоев по формуле
(5).
Описанная выше схема
аппроксимирует задачу с точностью до О( t +h2).
Невысокий порядок аппроксимации по 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 временных слоев с печатью результатов на
каждом слое. Таким образом, задача имеет вид
( ¶ 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. Программа, и результаты вычисления
приведены далее.
{
************************************************************* }
{ Приложение 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.