Модели на основе ОДУ. Краевая задача.

ЛЕКЦИЯ №7 МОДЕЛИ НА ОСНОВЕ ОДУ. КРАЕВАЯ ЗАДАЧА 1. Получение разностной схемы для линейного уравнения 2-го порядка с краевыми условиями 3-го рода интегро - интерполяционным методом Выше при построении разностной схемы нами применялся простой метод разностной аппроксимации, когда производные в уравнении и краевых условиях напрямую заменялись их разностными аналогами. В случае квазилинейных уравнений или в задачах с разрывными коэффициентами данный метод приводит к нарушению законов сохранения на сетке и появлению фиктивных источниковых слагаемых в разностном уравнении. Чтобы избежать появления указанных эффектов, применяют так называемый интегро - интерполяционный метод получения разностной схемы. Суть метода рассмотрим на примере решения уравнения d  du   k ( x)   p( x) u  f ( x)  0 dx  dx  (1) с краевыми условиями достаточно общего вида : слева - II рода, справа - III рода x  0,  k (0) du  F0 , dx , x  l ,  k (l ) du   (u (l )   ) dx где  ,  - известные числа. Введем сетку в области интегрирования уравнения [0, l]: h  xn : xn  nh, n  0,1,....N , h  l / N. Для построения разностной схемы выберем на сетке шаблон {xn-1, xn, xn+1} и ячейку {xn-1/2, xn+1/2} (рис.1). xn-1 xn xn+1 xn- 1/2 xn+ 1/2 Рис. 1. Шаблон и ячейка (затенена) на сетке Обозначим F  k ( x) du dx (2) По смыслу (2) - это поток. Интегрируем уравнение (1) с учетом (2) на ячейке xn 1 / 2 x x n 1 / 2 n 1 / 2 dF   dx   p( x) u dx   f ( x) dx  0 . dx xn 1 / 2 xn 1 / 2 xn 1 / 2 Выполняя точное интегрирование в первом слагаемом и применяя метод средних для приближенного численного вычисления остальных интегралов, получим  (Fn 1/ 2  Fn 1/ 2 )  pn ynh  fnh  0 , где pn  p( xn ), f n  f ( xn ) Проинтегрируем (2) на интервале [xn, xn+1] x n 1  xn du dx   dx x n 1  xn F dx . k ( x) Применяя метод средних справа, найдем (3) yn 1  yn   Fn 1 / 2 x n 1  xn dx , k ( x) или Fn1 / 2   n1 / 2 yn  yn1 , где  n1 / 2  h xn 1  xn h . dx k ( x) (4) Аналогично Fn 1 / 2   n 1 / 2 yn 1  yn ,  n 1 / 2  h h xn (5) dx x k ( x) n 1 Для величин  n  1/ 2 можно получить различные приближенные выражения, численно вычисляя интеграл методом трапеций или методом средних. Имеем, соответственно, две формулы, дающие близкие результаты  n 1 / 2  2kn kn 1 k k ,  n 1 / 2  kn 1 / 2  n n 1 . 2 k n  k n 1 (6) Аналогично  n 1 / 2  2k n k n 1 k k ,  n 1 / 2  kn 1 / 2  n n 1 . 2 k n  k n 1 Теперь, подставляя в (3) выражения для потоков (4), (5) и приводя подобные члены, получим окончательно систему линейных уравнений с трехдиагональной матрицей An yn1  Bn yn  Cn yn1   Dn , 1  n  N  1 , где An   n1/ 2 h , (7) Cn   n1/ 2 h , Bn  An  Cn  pnh , Dn  f n h . Система (7) совместно с краевыми условиями решается методом прогонки. Покажем далее, как тем же методом получают разностные аналоги краевых условий на примере краевого условия при x=0. Проинтегрируем уравнение (1) с учетом (2) на отрезке [0, x1|2]  x1 / 2  dF dx  dx x1 / 2 x1 / 2  p( x) u dx   f ( x) dx  0 . Второй и третий интегралы вычислим методом трапеций h h  ( F1 / 2  F0 )  ( p1 / 2 y1 / 2  p0 y0 )  ( f1 / 2  f 0 )  0 . 4 4 Далее, полагая y1/ 2  n=0 F1 / 2  1 / 2 y0  y1 и подставляя выражение для F1 / 2 согласно (4) при 2 y0  y1 , где 1 / 2  h x1 / 2  h , придем к формуле dx k ( x) M 0 y0  K 0 y1  P0 , (8) где h2 h2 M 0  1/ 2  p1/ 2  p0 , 8 4 . h2 h2 K0  p1/ 2  1/ 2 , P0  hF0  ( f1/ 2  f 0 ) 8 4 Для расчета величин в половинном узле можно принять простую аппроксимацию p1/ 2  p0  p1 , 2 f1/ 2  f 0  f1 . 2 Разностный аналог краевого условия при x  l получается аналогичным образом, если проинтегрировать уравнение (1) с учетом (2) на отрезке [xN-1/2, xN] и учесть, что поток FN  ( yN   ) , а FN 1 / 2   N 1 / 2 y N 1  y N . h В результате разностное краевое условие при x  l приводится к стандартному виду K N yN 1  M N yN  PN . (9) Сведем воедино полученные выше уравнения (7)-(9), составляющие разностную схему, аппроксимирующую исходную дифференциальную задачу (1) An yn1  Bn yn  Cn yn1   Dn , 1  n  N 1 M 0 y0  K 0 y1  P0 , K N yN 1  M N yN  PN . Отметим, что при уменьшении шага (в пределе при h  0 ), в (8) членами, содержащими h2 , можно пренебречь и тогда (8) преобразуется к виду y0  y1  т.е.  1 / 2 hF0 1 / 2 , y1  y0  F0 , что совпадает с выражением, которое можно получить, выh полняя простейшую аппроксимацию производной односторонней разностью, если брать χ в точке x1/2. В этом случае аппроксимация дает точность порядка O(h), тогда как (8) обеспечивает точность O(h2), совпадающую с порядком точности системы (7). Порядок точности, совпадающий с порядком аппроксимации, вообще говоря, надо проверять отдельно, оценивая невязку для разностного аналога краевого условия путем раз- ложения в ряды Тейлора. Процедура такого оценивания будет приведена позже при изучении уравнений в частных производных. В заключение данного раздела очертим алгоритм решения выписанной выше разностной схемы. Будем использовать для иллюстрации общих подходов простейшую аппроксимацию первых производных в краевых условиях односторонними разностями с порядком точности O(h). (В реальном моделировании строить разностные аналоги краевых условий надо так, как это сделано выше). В итоге, на левом краю при x=0 получим  k0 y1  y0  F0 . h (10) На правом краю при x=l запишем  kN y N  y N 1  ( yN   ) . h (11) Из (10) y0  y1  hF0 , k0 и, вспоминая основную прогоночную формулу yn   n1 yn1  n1 , записанную при n=0, т.е. y0  1 y1 1 , найдем начальные значения прогоночных коэффициентов 1 1, 1  hF0 . k0 Далее по известным рекуррентным формулам определим все прогоночные коэффициенты до последнего узла n=N. Наконец, подставляя в уравнение (11) выражение y N 1   N y N  N , получим yN   N k N  h . k N (1   N )  h (12) Теперь в обратном ходе окончательно находим по основной прогоночной формуле искомые значения сеточной функции y n В начале лекции при формулировке краевой задачи на границе x  l было поставлено краевое условие III рода. Это условие имеет достаточно общий характер, в частности, при k N  0,   1 оно переходит в краевое условие I рода - u   . Видно, что при k N  0,   1 формула (12) также переходит в выражение y N   , как и должно быть в соответствующей разностной схеме. .