МУ к лабораторным работам (1238837), страница 20
Текст из файла (страница 20)
Поэтому логично дляполучения одного недостающего уравнения аппроксимировать сжелаемым порядком точности дифференциальное уравнение дляq из (10.7) в ячейках сетки, примыкающих к левой границе.Аналогичным образом второе дополнительное соотношение получается аппроксимацией уравнения для r из (10.7) по ячейке,примыкающей к правой границе.За м еча ние. В схеме СХИ1 разностные уравнения для q приm = 0 и для r при m = M, можно рассматривать как аппроксимации первого порядка характеристических соотношений: для qвблизи левой границы и для r вблизи правой.В этом случае разностные уравнения (10.11) аппроксимируют соответствующие дифференциальные уравнения для инва144риантов со вторым порядком точности. Поэтому желательно,чтобы необходимые замыкающие соотношения обеспечивали туже точность.Приведемсоответствующий пример явной аппроксимаBции дифференциального уравAнения для q вблизи левой граниCцы.0xНарис.
10 изображенфрагментсеточной областиРис. 10вблизи левой границы ( x 0).Точка А является точкой пересечения c — характеристики,выпущенной из узла B с координатами (t p 1 , x0 0), с предыдущим слоем по времени.Интегрируя характеристическое соотношениеdqdt cFпо отрезку AB, получимqB q A F dt.ABЗаменим интеграл, например, по формуле трапеций. Имеем:qB q A 2( F A FB ) O( 3 ).Значение q A определим по известным значениям q в узлах p-гослоя по времени с помощью интерполяционной формулы:ppq A q0 ppppq1 q 0q 2q1 q 0 1 ( h ) 2 O( 3 ),2hh2здесь = a — расстояние точки A от левой границы.Подставляя выражение для q A в предыдущую формулу, получаем одно из недостающих дополнительных соотношений:p 1q0p q B q0 pppppq1 q 0q 2 q1 q 0 1 ( h ) 22hh2145 [ F (t p1 , 0) F (t p , a)] O( 3 3 ).2Другое недостающее уравнение получим, аппроксимируяаналогичным образом характеристическое соотношение для инварианта r вблизи правой границы:p 1rMppppprM rM 1r 2rM 1 rM 2 1 ( h) M2hh2p rM [ F (t p1 , 1) F (t p , 1 a)] O( 3 3 ).2Расчет по схеме СХИ2.
Из соотношений (10.13) находятсяq и r во всех точках начального слоя. Затем в рамках стандартного программного цикла (для p = 0, 1, 2, …, P – 1) осуществляется расчет данных на ( p 1) -м слое по времени. При этом:1) из уравнений (10.11) находится q mp 1 , а также rmp 1для m = 1, 2, …, M – 1;2) из уравнений (10.12) с использованием найденных издополнительных соотношений q0p 1 и rMp 1 отыскиp 1p 1и qM.Теперь рассмотрим примеры явных разностных схем непосредственно для задачи (10.6).ваются r0Схема СХ1. В этой схеме используется шаблонp + 1, m**p, m 1**p, mp, m + 1Разностные уравнения для внутренних узлов сетки:p 1ump 1vmapv m146ppumapv m 1 v m 12hp Fm pppa 2 u m1 2u m u m 1,2h2pppppu m1 u m 1 a 2 v m 1 2v m v m1,2h2h2p = 1, 2, …, P – 1;m = 1, 2, …, M – 1.(10.14)Аппроксимация краевых условий:u (t , 0) 1(t ),u(t , 1) 2 (t ),(10.15)p = 1, 2, …, P – 1.Начальные данные:0 ( x ),um1 m0 0,vmm = 0, 1, …, M.(10.16)Как и в случае схемы СХИ2 данная система соотношенийпока не замкнута.
Подходящий способ замыкания данной схемыобсуждается ниже. При разумном способе замыкания схемаимеет первый порядок точности и устойчива, если выполненоусловие Куранта: Cu a / h 1.Дополнительные соотношения для схемы СХ1. Недостающие для вычисления значений в точках левой и правой границысоотношения можно получить, аппроксимируя, например, дифференциальное уравнение для v по тому или иному шаблонувблизи левой и правой границ. Возможны следующие варианты.A) При m = 0: шаблон «явный правый уголок», уравнениеp 1pv0v 0ppu1 u 0 0.haПри m = M: шаблон «явный левый уголок», уравнениеp 1pvMv Mppu M u M 1 0.haДополнительные соотношения:p 1v0pv0 appи(u u 0 )h 1p 1vMpvM ahpp(u M u M 1 ).B) При m = 0: шаблон «неявный левый уголок», уравнениеp 1v0pv 0p 1au1p 1 u0h 0.При m = M: шаблон «явный правый уголок», уравнение147p 1pvMp 1v MauMp 1 u M 1h 0.Дополнительные соотношения:p 1v0p 1vMapv0 phvM p 1(u1ahp 1(u Mp 1 u0),p 1 u M 1 ).За м еча ние 1.
Возникают естественные вопросы: почему выбрано дифференциальное уравнение для v при конструировании дополнительных соотношений? Почему выбраны те илииные шаблоны? Точного ответа на первый вопрос нет. Что касается шаблона, то тут руководствуемся тем, чтобы не ухудшитьаппроксимацию и устойчивость схемы в целом. Вариант A), повидимому, непригоден, так как аппроксимация уравнений (10.7)по таким шаблонам, вообще говоря, неустойчива.За м еча ние 2. Данная схема является «близкой родственницей»схемы СХИ1 для инвариантов. В самом деле, уравнения группы(10.14) этой схемы являются линейной комбинацией (суммой иразностью) соответствующих уравнений схемы СХИ1.
Тем самым проясняется вопрос о происхождении данной схемы.Учитывая отмеченное обстоятельство, можно сделать вывод,что естественным способом замыкания данной схемы являетсяпривлечение в качестве дополнительных соотношений уравненийдля q (при m = 0) и r (при m = M) из группы уравнений (10.14)схемы СХИ1. В этом случае, очевидно, рассматриваемая схемаалгебраически эквивалентна схеме СХИ1 для инвариантов.Схема СХ2. В этой схеме используется шаблонp + 1, m**p, m 1148*p, m*p, m + 1Разностные уравнения для внутренних узлов сетки представляют собой линейные комбинации (полусумму и полуразность) уравнений для q и r схемы СХИ2 для инвариантов:p 1umpp umapv m 1 v m12hp Fm p 1vmpppppu 2u m u m1 p( Ft) m, a 2 m 12h2puua a m1 m 1 2h2v mp = 1, 2, …, P – 1;pppv m 1 2v m v m1 p( Fx ) m,ah2m = 1, 2, …, M – 1.(10.17)Краевые условия (10.15) и начальные данные (10.16) записываются так же как в схеме СХ1.Как и в случае схемы СХИ2 для инвариантов, данная система соотношений пока не замкнута.
Подходящий способ замыкания обсуждается ниже. При соответствующем способе замыкания схема имеет второй порядок точности и устойчива, есливыполнено условие Куранта: Cu a / h 1.Дополнительные соотношения для схемы СХ2. Недостающие для вычисления значений v в точках левой и правой границы соотношения можно получить, как и для схемы СХ1, аппроксимируя, например, дифференциальное уравнение для v потому или иному шаблону вблизи левой и правой границ.Способы, использованные для замыкания схемы СХ1, приводят к погрешности аппроксимации первого порядка, в то время, как уравнения для внутренних узлов данной схемы обеспечивают второй порядок аппроксимации.В качестве подходящего способа получения дополнительных уравнений, не портящих второй порядок точности схемы вцелом, применим аппроксимацию уравнения для v системы (10.6)вблизи левой и правой границы по прямоугольному шаблону.При m = 0 имеемp 1v0pv 0p 1v1p v1p 1au1p 1 u0hpapu1 u 0 0.h149Отсюда получаем дополнительную расчетную формулу, из которой находится v 0p 1.Аналогичным образом, при m = M получаем еще одну неp 1достающую формулу (для v M ).Применительно к данной схеме справедливы замечания,сделанные при обсуждении способов замыкания схемы СХ1.
Естественным способом представляется привлечение характеристических соотношений. В частности, если использовать в качестведополнительных уравнений аппроксимацию характеристическихсоотношений вблизи границ со вторым порядком точности, тоданная схема будет алгебраически эквивалентна СХ1.Расчет по схемам СХ1 и СХ2. Из соотношений (10.16) находятся u и v во всех точках начального слоя. Затем осуществляется расчет данных на (p + 1)-м слое по времени для значенийp = 0, …, P – 1.
При этом:1) из уравнений (10.14) находятся величины u mp 1 и v mp 1для m = 1, 2, …, M – 1;2) из (10.15) с использованием дополнительных соотношеp 1ний отыскиваются u0p 1, v 0p 1p 1и uM, vM .10.10. Контрольные вопросы1. Исследуйте на аппроксимацию и устойчивость разностныесхемы, рассматриваемые в данной лабораторной работе:1) схему «крест»;2) неявную схему с шаблоном*****3) неявную схему с шаблоном*****150**4) явную схему с шаблоном*****5) схему СХИ1 для инвариантов;6) схему СХИ2 для инвариантов;7) схему СХ1 для системы (10.6);8) схему СХ2 для системы (10.6).10.11. Порядок выполнения работы1. Проследите, как воспроизводится по различным разностнымсхемам решение «эталонных» задач (с известным точным решением).За м еча ние. Точное решение можно либо получить в аналитической форме самостоятельно, либо можно посмотреть на графическое его воспроизведение, выбрав нужную задачу по пунктам меню: «Методы/Аналитическое решение», и осуществив затем пуск соответствующей программы визуализации из пакетаменю «Запуск».2.
Объясните поведение численного решения «эталонной задачи» номер 2, полученного по разным схемам.3. Подберите параметры исходной задачи так, чтобы решениепредставляло собой стоячую волну. Проследите, как это решение воспроизводится по различным разностным схемам.Под стоячей волной здесь подразумевается стационарноерешение волнового уравнения. Например, для задачи 2 sin x,u tt u xxu ( 0, t ) u (1, t ) 0,u ( x, 0) sin xможно рассмотреть решение u ( x, t ) sin x.10.12.
Библиографическая справкаО простейших схемах для решения волнового уравнения можно прочитать в [1–4, 27, 36]. О сеточно-характеристическихметодах см. [35].151Л АБОР АТ ОРН АЯ РАБОТ А 11ЧИСЛЕННОЕ РЕШЕНИЕДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙВ ЧАСТНЫХ ПРОИЗВОДНЫХПАРАБОЛИЧЕСКОГО ТИПА.УРАВНЕНИЕ ТЕПЛОПРОВОДНОСТИ11.1. ВведениеЭта работа знакомит с численными методами решения дифференциальных уравнений в частных производных параболического типа на примере одномерного линейного уравнения теплопроводности.
Рассматривается следующая краевая задача: f (t , x),u t a 2 u xxt 0,u (0, x) ( x ),x1 x x 2 ,u (t , x1 ) 1 (t ),t 0,u (t , x2 ) 2 (t ),t 0.x1 x x 2 ,На примерах рассматриваются особенности поведения численных решений этого уравнения в зависимости от входныхданных задачи, выбора разностного метода, сеточных параметров (шагов по времени и координате).Реализованы следующие разностные схемы для численногорешения одномерного уравнения теплопроводности:1) шеститочечная параметрическая схема;2) схема Франкела–Дюфорта;3) схема Ричардсона;4) явная центральная четырехточечная схема;5) схема Алена–Чена;6) нецентральная явная четырехточечная схема;7) схема Саульева.Полученные приближенные (численные) решения однородного уравнения можно сравнить с точными для нескольких тестовых задач.11.2.
Дифференциальная краевая задачаКак уже было отмечено, в работе рассматривается задача: f (t , x),u t a 2 u xxt 0,u (0, x) ( x ),x1 x x 2 ,u (t , x1 ) 1 (t ),t 0,u (t , x2 ) 2 (t ),t 0.x1 x x 2 ,11.3. Сеточная областьДля рассмотренной задачиW h [(t p , xm )],pu h [u m ],p = 0, 1, ..., P,p = 0, 1, ..., P,m = 0, 1, …, M,m = 0, 1, …, M,где u mp — компонента сеточной функции, относящаяся к узлу(t p , x m ), t p p, — шаг по времени, P T , x m x1 mh,h — шаг по координате, Mh x2 x1.11.4. Пример разностной задачи (разностной схемы)Для рассмотренной дифференциальной задачи одна из возможных разностных схем имеет следующий вид:p 1um a2ppu m1 2u m u m1h2p fm ,p = 0, 1, …, P – 1;m = 1, 2, …, M – 1;0 ,ummm = 0, 1, …, M;ppu0 1 ,ppuM 2 ,152pp ump = 1, 2, …, P;p = 1, 2, …, P.(11.1a)11.5.