МУ к лабораторным работам (1238837), страница 19
Текст из файла (страница 19)
Разностные уравнения для внутренних узлов сетки:p 1ump2pp 1 2u m u m a2ppu m 1 2u m u m1p fm ,2h(10.2)p = 1, 2, …, P – 1;m = 1, 2, …, M – 1.Аппроксимация краевых условий:pu0 1(t p ),u mp 2 (t p ),p = 1, 2, …, P.Аппроксимация начальных условий:0 ( x ),um1 m0u 1m u m 2 ( x m ),m = 0, 1, …, M.Во внутренних узлах уравнения (10.2) аппроксимируют исходp + 1, mное дифференциальное уравнение со*вторым порядком точности. Однако вданном варианте схемы второе на**чальное условие аппроксимируется p*,m 1p, mp, m + 1простейшим образом — с первым поp *1, mрядком точности (по ).
Поэтому в целом это схема первого порядка.Условие устойчивости численного решения — число Куранта Cu a / h 1. (По поводу отмеченных фактов см. п. 10.4.–10.7). О проведении расчетов по этой схеме см. п. 10.8.136Явная схема. Шаблон схемы имеет вид:*****Разностные уравнения для внутренних узлов сетки:p 1umpp 1p 1 2u m u m a22p 1u m 1 2u mp 1 u m1h2p 1 fm,(10.3)p = 1, 2, …, P – 1;m = 1, 2, …, M – 1.Начальные и краевые условия аппроксимируются как и в схеме«крест». Таким образом, данная схема состоит из уравнений(10.3) и начальных и краевых условий из схемы «крест».Это схема первого порядка точности, но абсолютно неустойчива! Для решения конкретных задач она не используетсяи приводится здесь лишь как пример абсолютно неустойчивойразностной схемы.Неявная схема (1).
Шаблон схемы имеет вид:*****Разностные уравнения для внутренних узлов сетки:p 1umpp 1p 1 2u m u m2 a2p 1u m1 2u mp 1 u m 1h2p 1 fm,(10.4)p = 1, 2, …, P – 1;m = 1, 2, …, M – 1.Начальные и краевые условия аппроксимируются также как всхеме «крест». Таким образом, данная схема состоит из приведенных здесь уравнений (10.4) и начальных и краевых условийиз схемы «крест».Это схема первого порядка точности абсолютно устойчива.
Однако для решения конкретных задач она практически не137используется, так как дает слишком плохие результаты. В этомможно убедиться на методических примерах этой лабораторной работы.Неявная схема (2). Шаблон схемы имеет вид:*******Разностные уравнения для внутренних узлов сетки:p 1umpp 1 2u m u m2p 1 a2p 1(u m1 2u mp 1p 1p 1 u m 1 ) (u m 1 2u m2h 2p 1 u m 1 )p fm ,(10.5)p = 1, 2, …, P – 1;m = 1, 2, …, M – 1.Начальные и краевые условия аппроксимируются как в схеме«крест». Таким образом, данная схема состоит из уравнений(10.5) и начальных и краевых условий из схемы «крест».Это схема первого порядка точности по t (за счет грубойаппроксимации начальных данных, как и схема «крест»), абсолютно устойчива.10.8.2.
Об аппроксимации начальных данных. В поясненияхк схеме «крест» приведена простейшая аппроксимация условияut (0, x) 2 (t ), приводящая к погрешности O() , в то времякак во внутренних узлах сетки погрешность аппроксимации длянекоторых из рассматриваемых здесь схем является величинойвторого порядка как по , так и по h.Можно без труда обеспечить второй порядок и при аппроксимации данного начального условия, так что соответствующаясхема станет схемой второго порядка точности.Представим значение решения в точках первого слоя повремени в виде ряда Тейлора по степеням :2u (, x m ) u (0, x m ) u t (0, x m ) u tt (0, x m ) O( 3 ).2138Замечаем, что из дифференциального уравнения следует f (t , x ).u tt a 2 u xxТаким образом,u (, x m ) u (0, x m ) u t (0, x m ) 2 (0, x m )} O(3 ). { f (0, xm ) a 2u xx2Отсюда получаем расчетную формулу для данных на первомслое по времени:2u 1m 1 ( x m ) 2 ( x m ) { a 2 1 xx ( x m ) f (0, x m )}.2Последнее соотношение, переписанное в виде0u 1m u m 2 ( x m ) {a 2 1 xx ( x m ) f (0, x m )},2очевидно, аппроксимирует условие ut (0, x ) 2 (t ) со вторымпорядком точности.10.8.3.
О последовательности вычислений. Сначала из соотношений для начальных данных схемы «крест» находятся значе0 (m = 0, …, M) и u1 (m = 1, 2, …, M – 1). Недостающиения: u mmкомпоненты сеточной функции на первом слое u10 и u1M определяются из условий на границах. Затем для p = 1, 2, …, P – 1 изсоотношения (10.2) отыскиваются u mp 1 (m = 0, 1, …, M).
Длянеявных схем необходимо решать систему линейных уравненийс трехдиагональной матрицей.10.9. Сведение задачи (10.1) к задачедля системы двух уравнений первого порядкаЗа м еча ние. В этом пункте демонстрируется, как с помощьюискусственного приема рассматриваемую задачу можно свести кдругой — известной как задача для уравнений акустики. Последние — ничто иное, как уравнения переноса. Таким образом,этот раздел представляет собой связующее звено между даннойлабораторной работой и работой «Численное решение диффе-139ренциальных уравнений в частных производных гиперболического типа. Уравнение переноса».Задача (10.1) эквивалентна задаче:ut av x F (t , x),v t au x 0,t 0,0 x 1,где функция, стоящая в правой части первого уравненияtF (t , x ) f (, x) d 2 ( x ),0u (0, x ) 1( x ),v (0, x) 0,u (t , 0) 1(t ),u (t , 1) 2 (t ).(10.6)10.9.1.
Обоснование. Введем в рассмотрение новую функциюv (t , x) связав ее с u (t , x) уравнением v t au x 0 и полагаяv (0, x) 0. Тогда волновое уравнение можно записать в виде f (t , x) 0.utt av txИнтегрируем последнее соотношение по t, получаем:t [(ut av x )t f (t , x)] dt 0.0После выполнения интегрированияtu t av x f (t , x) dt [u t av x ] t 0 ,0tт. е. ut av x F (t , x), где F (t , x) f (, x) d 2 ( x).0Таким образом, задача (10.1) сведена к задаче (10.6) длядвух дифференциальных уравнений первого порядка.10.9.2.
Дополнительные замечания. Легко проверить непосредственно, что дифференциальные уравнения (10.6) можнозаписать в виде:qt aqx F (t , x),rt a rx F (t , x ),где q u v , r u v — инварианты Римана.140Особенность последних уравнений состоит в том, что в каждом из них дифференцируется лишь одна неизвестная функция(q или r). Отметим, что возможность записи исходной системытипа (10.6) в виде системы для инвариантов (функций, для которых уравнения разделяются так же, как и здесь) является признаком гиперболичности системы.Очевидно, что рассматриваемая задача может быть сформулирована как задача для инвариантов:qt aqx F (t , x),q(0, x ) 1 ( x),rt arx F (t , x ),t 0,r (0, x) 1( x ),q(t , 0) r (t , 0) 21 (t ),0 x 1,(10.7)q(t , 1) r (t , 1) 2 2 (t ).Уравнения для инвариантов, к которым мы перешли, являются известными уравнениями переноса (см.
лабораторную работу «Численное решение дифференциальных уравнений в частных производных гиперболического типа. Уравнение переноса»). Последнее обстоятельство позволит распространитьразностные схемы для уравнения переноса на наш случай.Таким образом, вместо исходной задачи (10.1) можно решать задачу (10.6) или задачу (10.7) для инвариантов. В последнем случае исходная искомая функция u вычисляется через q и rв нужных точках по формуле u (q r ) / 2.10.9.3.
Варианты разностных схем для задачи (10.6). Предварительные замечания. Как известно, неявные схемы практически всегда устойчивы. Однако, применительно к задаче (10.6)использование их связано с определенными трудностями, таккак расчет данных на очередном слое по времени требует решения линейной системы уравнений (более общего вида, нежелисистема с трехдиагональной матрицей). Поэтому в рамках данной работы мы сосредоточим внимание на способах построенияболее простых (с точки зрения реализации) явных разностныхсхем для задачи (10.6). С ними, однако, применительно к задаче(10.6) снова не все ясно. (Попробуйте по произвольному «явному» шаблону построить устойчивую схему!).
Поэтому сначаларазберемся как строятся подходящие явные схемы для задачи(10.7), сформулированной для инвариантов q и r.Так как каждое из дифференциальных уравнений (10.7) естьуравнение переноса, то соответствующие явные схемы для последнего с некоторыми уточнениями, касающимися расчета в точкахправой и левой границы, легко обобщаются на задачу (10.7).141СХИ1 — схема для инвариантов (первого порядка точности). Эта схема является естественным обобщением на случайсистемы для инвариантов схем типа «явный угол» для уравнения переноса.
С учетом характеристических направлений (направлений переноса инвариантов) первое из уравнений (10.9)аппроксимируется по шаблонуp + 1, m***p, mp, m + 1шаблон для второго уравненияp + 1, m**p, m 1*p, mРазностные уравнения для внутренних узлов сетки:p 1qmpp qmpq m 1 q mahp = 1, 2, …, P – 1;p 1rmpp rmap Fm ,m = 0, 1, …, M – 1;prm rm 1p Fm ,h(10.8)p = 1, 2, …, P – 1;m = 1, 2, …, M.Аппроксимация краевых условий:q(t , 0) r (t , 0) 21 (t ),(10.9)q(t , 1) r (t , 1) 2 2 (t ),p = 1, 2, …, P – 1.Начальные данные:0 ( x ),qm1 mrm0 1 ( x m ),m = 0, 1, …, M.(10.10)Схема имеет первый порядок точности, устойчива при выполнении условия Куранта Cu a / h 1.Последовательность вычислений для схемы СХИ1.
Из соотношений (10.10) находятся q и r во всех точках начальногослоя. Затем в рамках стандартного программного цикла для142p = 0, 1, 2, …, P – 1 осуществляется расчет данных на (p + 1)-мслое по времени. При этом:p 11) из соотношений (10.8) находятся qmдля m = 0, …, M – 1и rmp 1 для m = 1, 2, …, M;2) из (10.9) с использованием найденных значений q0p 1 иp 1p 1отыскиваются r0rMp 1и qM .СХИ2 — схема для инвариантов (второго порядка точности). В этой схеме используется шаблонp + 1, m****p, m 1p, mp, m + 1За м еча ние. Данная схема является обобщением на случай системы (10.7) явной четырехточечной схемы второго порядка дляуравнения переноса.Разностные уравнения для внутренних узлов сетки:p 1qmpp qmapq m 1 q m12hp Fm p 1rmpp rmapppq 2 q m q m1 pp( Ft) m, a ( Fx ) m a 2 m12h2prm 1 rm12hr 2 rm rm 1 p pp, Fm ( Ft) m a ( Fx ) m a 2 m 12h2pp = 1, 2, …, P – 1;ppm = 1, 2, …, M – 1.(10.11)Аппроксимация краевых условий:q(t , 0) r (t , 0) 21 (t ),(10.12)q(t , 1) r (t , 1) 2 2 (t ),p = 1, 2, …, P – 1.143Начальные данные:0 ( x ),qm1 mrm0 1 ( xm ),m = 0, 1, …, M.(10.13)В отличие от схемы СХИ1 данная система соотношенийпока не замкнута.
Подходящий способ замыкания данной схемырассматривается ниже. При выбранном там способе замыканиясхема имеет второй порядок точности и устойчива при выполнении условия Куранта: Cu a / h 1.Дополнительные соотношения схемы СХИ2. Заметим, чтодифференциальное уравнение для инварианта q — ничто иное,как характеристическое соотношениеdqdt c F,выполняющееся вдоль характеристик c : dx / dt a, которыепредставляют собой семейство прямых x at const.В частности, это характеристическое соотношение связывает значение q в точках левой границы со значениями q внутрирасчетной области (на c — характеристике, выпущенной израссматриваемой точки левой границы).