Реферат: Решение смешанной задачи для уравнения гиперболического типа методом сеток
Название: Решение смешанной задачи для уравнения гиперболического типа методом сеток Раздел: Рефераты по математике Тип: реферат | |||||
Решение смешанной задачи для уравнения гиперболического типа методом сеток Рассмотрим смешанную задачу для волнового уравнения ( ¶ 2 u/ ¶ t2 ) = c 2 * ( ¶ 2 u/ ¶ 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 ) = ( ¶ 2 u/ ¶ x2 ), то в дальнейшем будем считать с = 1. Для построения разностной схемы решения задачи строим в области в = {(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) .
(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. Программа, и результаты вычисления приведены далее. Приложение 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 {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. Приложение 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 {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. (выполнения лабораторной работы. Вариант 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 {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. ( выполнения лабораторной работы. Вариант 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 {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. ( выполнения лабораторной работы. Вариант 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 {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. ( выполнения лабораторной работы. Вариант 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 {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. |