Определить порядок аппроксимации схемы лакса. Примеры разностных схем. Явные и неявные разностные схемы. Порядок точности разностной схемы

Различают три метода построения разностных схем на заданном шаблоне:

· метод разностной аппроксимации;

· интегро-интерполяционный метод;

· метод неопределенных коэффициентов.

Методом разностной аппроксимации мы уже пользовались при составлении схем (24), (26). Согласно данному методу, каждая производная, входящая в уравнение и краевое условие, заменяется каким-либо разностным выражением с учетов узлов заданного шаблона. Метод позволяет легко составить разностные схемы с первым и вторым порядком аппроксимации, когда коэффициенты уравнения достаточно гладкие функции. Обобщение данного подхода на ряд важных случаев затруднителен. Например, если коэффициенты уравнения разрывные, или предполагается использовать непрямоугольную и неравномерную сетку, возникает неопределенность в построении разностной схемы.

При использовании интегро-интерполяционного метода или метода баланса используют дополнительные физические соображения, сводящиеся к составлению уравнений сохранения тех или иных величин. В данном методе после выбора шаблона область разбивается на ячейки. Дифференциальное уравнение интегрируют по ячейке и по формулам векторного анализа приводят к интегральной форме, отвечающей некоторому интегральному закону. Интегралы вычисляют приближенно по одной из квадратурных формул и получают разностную схему.

Представим уравнение теплопроводности с переменным коэффициентом теплопроводности в виде: . Выберем для его аппроксимации шаблон, представленный на рис.8, где пунктиром выделена соответствующая ячейка.

Выполним интегрирование по ячейке:

и аппроксимируем первый интеграл по формуле средних, а второй интеграл - по формуле прямоугольников, тогда

В последнем выражении производные заменим конечными разностями и, считая сетку равномерной, получим разностную схему

Если k = const, то схема (35) совпадает с неявной схемой (24).

Рис.8. Шаблон и ячейка интегро-интерполяционного
метода для уравнения теплопроводности

Интегро-интерполяционный метод наиболее полезен, когда коэффициенты уравнения является негладкими или даже разрывными. В этом случае обращение к более общим - интегральным законам возвращает нас к более правильным обобщенным решениям.

Рассмотрим пример использования разностной схемы (35) для расчета теплопроводности среды, состоящей из трех сред с разными коэффициентами теплопроводности, т.е.

(36)

где k 1 , k 2 , k 3 , вообще говоря, различные неотрицательные числа. Исходное уравнение можно в этом случае записать в виде:

(37)

Для расчета по схеме (35) с коэффициентом теплопроводности (36) будем полагать, что

а на левой x = 0 и правой x = a границе согласно (37) будем поддерживать нуль температуры, т.е. и .

На листинге_№4 приведен код программы, которая решает уравнение (36), (37) согласно разностной схеме (35), (38).

Листинг_№4

%Программа решения уравнения теплопроводности

%(37) с разрывным коэффициентом

%теплопроводности (36)

global a k1 k2 k3

%определяем отрезок интегрирования и

%три значения коэффициента теплопроводности

%в трех областях отрезка интегрирования

a=3; k1=0.1; k2=100; k3=10;

%определяем шаг по времени и по пространству

tau=0.05; h=0.05;

x=0:h:a; N=length(x);

%Строим начальное распределение температуры

if x(i)<=0.5*a

y(i)=((2*Tm)/a)*x(i);

if x(i)>0.5*a

y(i)=((2*Tm)/a)*(a-x(i));

%рисуем начальный профиль температуры

%толстой красной линией

plot(x,y,"Color","red","LineWidth",3);

%вычисляем коэффициенты прогонки A(n), B(n)

%C(n): A(n)y2(n+1)+B(n)y2(n)+C(n)y2(n-1)=y(n)

A(n)=-(tau/h^2)*k(x(n)+0.5*h);

B(n)=1+(tau/h^2)*...

(k(x(n)+0.5*h)+k(x(n)-0.5*h));

C(n)=-(tau/h^2)*k(x(n)-0.5*h);

%определяем левое граничное условие

alpha(2)=0; beta(2)=0;

alpha(n+1)=-A(n)/(B(n)+C(n)*alpha(n));

beta(n+1)=(y(n)-C(n)*beta(n))/...

(B(n)+C(n)*alpha(n));

%задаем правое граничное условие

for n=(N-1):-1:1

y(n)=alpha(n+1)*y(n+1)+beta(n+1);

%рисуем текущий профиль температуры

%определяем коэффициент теплопроводности

global a k1 k2 k3

if (x>=0)&(x<=a/3)

if (x>a/3)&(x<=(2*a)/3)

if (x>(2*a)/3)&(x<=a)

На рис.9 приведен итог работы кода программы листинга_№4. Красной жирной линией нарисован начальный треугольный профиль температуры. Вертикальные стрелки на графике отделяют области с разными коэффициентами теплопроводности. Согласно коду листинга_№4, коэффициенты теплопроводности отличаются друг от друга на три порядка.

Рис.9. Решение уравнения теплопроводности (37) с разрывным
коэффициентом теплопроводности (36)

Метод неопределенных коэффициентов заключается в том, что в качестве разностной схемы берут линейную комбинацию решений в узлах некоторого шаблона. Коэффициенты линейной комбинации определяют из условия максимального порядка соответствующей невязки по t и h .

Так, для уравнения на шаблоне рис.8 можем записать следующую схему с неопределенными коэффициентами

Определяем невязку

Подставим (31) в (40), тогда

(41)

Большинство членов в (41) обнуляются при условии

. (42)

Подставляя (42) в (39) получим разностную схему (24).

Метод неопределенных коэффициентов применим и к более сложным случаям. Например, для треугольной сетки, шаблон которой приведен на рис.10 можно получить следующую разностную схему

Рис.10. Шаблон треугольной сетки для разностного уравнения (43)

Рассмотрим нерегулярные узлы разностной схемы, т.е. ее граничные условия. Для уравнения теплопроводности u t = k u xx нерегулярными являются граничные узлы n = 0 и n = N . Если рассматривается первая краевая задача

то легко записать соответствующие разностные условия

которые выполняются точно, т.к. невязка для них равна нулю.

Более сложным является случай второй краевой задачи, когда граничное условие содержит производную по x . Например, при задании на краях теплового потока граничные условия приобретают следующий вид:

Производные в (44) можно аппроксимировать правой (левой) конечной разностью:

Невязка разностных уравнений (45) легко оценивается:

(46)

Таким образом, согласно (46), невязка граничных условий имеет первый порядок точности по h , тогда как в регулярных точках порядок точности второй по h , т.е. при выборе аппроксимации граничных условий по формулам (45) происходит потеря точности.

Для повышения точности граничных условий рассмотрим метод фиктивных точек . Введем вне отрезка две фиктивные точки: , и запишем в точках n = 0 и n = N явную разностную схему (26), тогда

Аппроксимируем левое и правое граничное условие с помощью центральной разности, т.е.

Исключая из (47), (48) фиктивные точки и значения функции в них, находим граничные условия второго порядка точности по h :

(49)

Граничные условия (49) являются явными, т.к. содержат только по одному значению на следующем слое.

Помимо метода фиктивных точек есть другой метод уменьшения невязки, он более универсален, но менее нагляден. Разложим u (t ,x 1) в окрестности x 0 , тогда

Согласно (44), , а из уравнения теплопроводности найдем . Подставляя данные оценки в разложение Тейлора, находим

Делая в (50) замену , получим левое граничное условие (49).

Согласно приведенной выше процедуре можно добиться повышенной точности в аппроксимации граничных условий.

Аппроксимация

Пусть задана область G переменных x = (x 1 ,x 2 ,…,x p ) с границей G и поставлена корректная задача решения уравнения с граничными условиями:

Au (x ) - f (x ) = 0, x Î G ; (51)

Ru (x ) - m (x ) = 0, x Î G. (52)

Введем в области G + G сетку с шагом h , которая содержит регулярные (внутренние) узлы w h и нерегулярные (граничные) узлы g h .

Перейдем в (51), (52) к соответствующим разностным аналогам

A h y h (x ) - j h (x ) = 0, x Î w h ; (51¢)

R h y h (x ) - c h (x ) = 0, x Î g h . (52¢)

Близость разностной схемы (51¢), (52¢) к исходной задаче (51), (52) определяется величинами невязок:

Разностная схема (51¢), (52¢) аппроксимирует задачу (51), (52), когда

аппроксимация имеет p -й порядок, когда

Дадим некоторые комментарии к выбору норм. Для простоты будем рассматривать одномерный случай, т.е. G = [a ,b ].

Можно использовать чебышевскую или локальную норму

,

или гильбертову, среднеквадратическую:

.

Часто строят ассоциированные или связанные с оператором A энергетические нормы. Например,

Выбор нормы регулируется двумя противоположными соображениями. С одной стороны, желательно, чтобы разностное решение y было близко к точному решению в наиболее сильной© норме. Например, в задачах на разрушение конструкций малость деформаций в не гарантирует целостность конструкций, а малость в норме - гарантирует. С другой стороны, чем слабее норма , тем легче разностную схему построить и доказать ее сходимость.

Функции y h , j h , c h , входящие в (51¢), (52¢), определены на сетке, поэтому для них необходимо определить соответствующие сеточные нормы , и . Обычно их вводят так, чтобы они переходили в выбранные нормы , и при h ® 0. В качестве разностных аналогов чебышевской и гильбертовых норм выбирают выражения

или близкие аналоги.

Устойчивость

Под устойчивостью (неустойчивостью) разностной схемы понимается то, что малые ошибки, возникающие в процессе счета (или внесенные с входными данными), при последующих расчетах уменьшаются (возрастают).

Рассмотрим пример неустойчивой разностной схемы для задачи Коши дифференциального уравнения u ¢ = a u . Выберем следующее однопараметрическое семейство разностных схем:

. (53)

Исследуем рост ошибки dy n начальных данных уравнения (53). Поскольку уравнение (53) линейно, постольку ошибка dy n удовлетворяет тому же уравнению (53). Изучим специальный вид ошибки dy n = l n . Подставим это представление в (53), тогда

Решение квадратного уравнения (54) при h ® 0 дает следующие оценки корней

Из оценок корней в (55) следует, что при s < ½ второй корень |l 2 | > 1, т.е. за один шаг ошибка возрастает в несколько раз. Проверим это.

На листинге_№5 приведен код программы, иллюстрирующей расчет по неустойчивой при s = 0,25 схеме (53) и по устойчивой схеме при s = 0,75. В начальных данных выбирались малые возмущения. Далее проводились серии расчетов с уменьшающимся значением шага сетки h . На рис.11 приведены итоговые графики зависимости значения возмущения в начальных данных на правом конце отрезка интегрирования в зависимости от шага сетки. Отчетливо видно сколь разительно отличаются друг от друга расчеты по неустойчивой и устойчивой схемам. Используя данную программу можно убедиться в пороговом значении параметра s = 0,5: при s < 0,5 схема неустойчива, при s ³ 0,5 - устойчива.

Листинг_№5

%Программа расчета по неустойчивой схеме при

%sigma=0,25 и по устойчивой схеме при sigma=0,75

%очищаем рабочее пространство

%определяем константу уравнения u"=alpha*u

%определяем значения sigma=0,25; 0,75

sigm=0.25:0.5:0.75;

for s=1:length(sigm)

%определяем начальное значение шага сетки

x=0:h:1; N=length(x);

%определяем возмущения начальных данных

dy(1)=1e-6; dy(2)=1e-6;

%осуществляем расчет возмущения начальных

%данных на правом конце отрезка интегрирования

dy(n+1)=(2+(alpha*h-1)/sigma)*dy(n)+...

(1/sigma-1)*dy(n-1);

%запоминаем возмущение на правом конце и

%шаг сетки

deltay(i)=dy(N);

%рисуем график зависимости возмущения на

%правой границе от шага сетки

plot(step,deltay);

Рис.11. Графики зависимости возмущения при расчете по
схеме (53) на правой границе от шага сетки h

Разностная схема (51¢), (52¢) устойчива , если решение системы разностных уравнений непрерывно зависит от входных данных j , c и эта зависимость равномерна относительно шага сетки. Уточним непрерывную зависимость. Это означает, что для любого e > 0 найдется такое d (e ), не зависящее от h , что

, (56)

Если разностная схема (51¢), (52¢) линейна, то разностное решение линейно зависит от входных данных. В этом случае можно положить, что d (e ) = e /(M + M 1), где M , M 1 - некоторые неотрицательные величины, независящие от h . В итоге условие устойчивости для линейных разностных схем можно записать в виде:

Непрерывную зависимость разностного решения от j называют устойчивостью по правой части , а от c - устойчивостью по граничным данным .

В дальнейшем будем рассматривать двуслойные разностные схемы , т.е. такие схемы которые содержат один известный и один новый, неизвестный слой.

Двуслойная разностная схема называется равномерно устойчивой по начальным данным, если при выборе начальных данных с любого слоя t * (t 0 £ t * < T ) разностная схема устойчива по ним, причем устойчивость равномерна по t * . Для линейных схем условие равномерной устойчивости можно записать в виде

где константа K не зависит от t * и h , - решения разностной схемы A h y = j с начальными данными и с одной и той же правой частью.

Достаточный признак равномерной устойчивости. Для равномерной устойчивости по начальным данным достаточно, чтобы при всех m выполнялось

Доказательство. Условие (60) означает, что если на некотором слое возникла ошибка dy , то при переходе на следующий слой норма возмущения ||dy || возрастает не более чем в (1 + Сt ) £ e C t раз. Согласно (59), при переходе со слоя t * на слой t требуется m = (t - t *)/t шагов по времени, т.е. ошибка возрастает не более чем в . В итоге имеем

что, согласно определению в (59), означает равномерную устойчивость по начальным данным.

Теорема. Пусть двухслойная разностная схема A h y = j равномерно устойчива по начальным данным и такова, что если два разностных решения A h y k = j k равны на некотором слое, т.е. , то на следующем слое выполняется соотношение

где a = const. Тогда разностная схема устойчива по правой части.

Доказательство. Помимо решения y рассмотрим решение , соответствующее возмущенной правой части . В дальнейшем будем считать, что . Это можно предположить, т.к. исследуется устойчивость по правой части.

конфигурация узлов, значения сеточной функции в которых определяют вид разностных уравнений во внутренних (не приграничных) точках сетки. Как правило, на рисунках с изображениями шаблонов точки, участвующие в вычислении производных, соединяются линиями.

Схема Куранта - Изаксона - Риса (КИР), которую иногда также связывают с именем С.К. Годунова, получается при , . Ее порядок аппроксимации . Схема КИР условно устойчива, т.е. при выполнении условия Куранта . Приведем разностные уравнения для схемы Куранта - Изаксона - Риса во внутренних точках расчетной области:

Эти схемы, имеющие также название схемы с разностями против потока (в англоязычной литературе - upwind) могут быть записаны в виде

Их преимущество состоит в более точном учете области зависимости решения. Если ввести обозначения

то обе схемы можно записать в следующих формах:

(потоковая форма разностного уравнения);

(здесь явно выделен член со второй разностью, придающий устойчивость схеме);

(уравнение в конечных приращениях).

Рассмотрим также метод неопределенных коэффициентов для построения разностной схемы правый уголок первого порядка точности для уравнения переноса

Схему можно представить в виде

Схема Куранта - Изаксона - Риса тесно связана с численными методами характеристик . Дадим краткое описание идеи таких методов.

Две последние полученные схемы (при разных знаках скорости переноса) можно интерпретировать следующим образом. Построим характеристику , проходящую через узел (t n + 1 , x m ), значение в котором необходимо определить, и пересекающую слой t n в точке . Для определенности считаем, что скорость переноса c положительна.

Проведя линейную интерполяцию между узлами x m - 1 и x m на нижнем слое по времени, получим

Далее перенесем вдоль характеристики значение u n (x") без изменения на верхний слой t n + 1 , т.е. положим . Последнее значение естественно считать приближенным решением однородного уравнения переноса. В таком случае

или, переходя от числа Куранта снова к сеточным параметрам,

т.е. другим способом пришли к уже известной схеме "левый уголок", устойчивой при . При точка пересечения характеристики , выходящей из узла (t n + 1 , x m , с n - м слоем по времени расположена левее узла (t n , x m - 1 ). Таким образом, для отыскания решения используется уже не интерполяция , а экстраполяция, которая оказывается неустойчивой.

Неустойчивость схемы "правый уголок" при c > 0 также очевидна. Для доказательства этого можно использовать либо спектральный признак, либо условие Куранта, Фридрихса и Леви. Аналогичные рассуждения можно провести и для случая c < 0 и схемы "правый уголок".


Неустойчивая четырехточечная схема получается при , ее порядок аппроксимации . Сеточные уравнения для разностной схемы будут иметь следующий вид:

Схема Лакса - Вендроффа возникает при . Порядок аппроксимации схемы Лакса - Вендроффа есть . Схема устойчива при выполнении условия Куранта .

Эту схему можно получить либо методом неопределенных коэффициентов, либо путем более точного учета главного члена погрешности аппроксимации. Рассмотрим процесс вывода схемы Лакса - Вендроффа подробнее. Проводя исследование предыдущей четырехточечной схемы на аппроксимацию (а исследование это довольно элементарно и сводится к разложению функции проекции на сетку точного решения дифференциальной задачи в ряд Тейлора), получим для главного члена погрешности

При выводе выражения для главного члена погрешности аппроксимации использовано следствие исходного дифференциального уравнения переноса

Которое получается путем дифференцирования исходного уравнения (3.3) сначала по времени t , затем по координате x и вычитанием одно из другого получившихся соотношений.

Далее, заменяя вторую производную во втором слагаемом в правой части с точностью до O(h 2) , получим новую разностную схему, аппроксимирующую исходное дифференциальное уравнение с точностью . Сеточные уравнения для схемы Лакса - Вендроффа во внутренних узлах расчетных сеток есть

Неявная шеститочечная схема возникает при q = 0 ; при ее порядок аппроксимации , при .

Пример 1. Разностная схема для уравнения Пуассона, относящегося к эллиптическому типу.

Рассмотрим построение разностной схемы для первой краевой задачи для уравнения А и = f(x,y) в области, представляющей собой прямоугольник со сторонами, параллельными осям координат. Пусть с этим прямоугольником связана равномерная сетка с шагами h x и h y .

Краевую задачу

можно записать в операторной форме:


Отметим, что в эту запись включены и граничные условия.

Заменив дифференциальные операторы разностными, получим уравнения


которые аппроксимируют исходное дифференциальное уравнение со вторым порядком 0(h 2 + h 2) точности и действуют во всех внутренних точках области.

Разностные аналоги краевых условий будут иметь вид

Разностная аппроксимация дифференциального уравнения совместно с разностными аналогами краевых условий образуют разностную схему для уравнения Пуассона.

Разностную схему можно по аналогии с краевой задачей записать в операторном виде:

где в L/, включены как разностное уравнение, так и разностное краевое условие:


Разностное уравнение связывает значения сеточной функции в пяти точках, образующих разностный шаблон для этого уравнения. Для данного случая этот шаблон называют крест. Можно представить и другие шаблоны для этого уравнения.

Мы получим приближенное решение дифференциальной краевой задачи, если определим значения сеточной функции во всех внутренних узлах области. Для этого необходимо решить совместно систему алгебраических линейных уравнений, размерность которой равна количеству внутренних узлов области. В этом случае говорят о неявной разностной схеме. Любое интересующее нас значение Uij можно определить лишь из решения всей разностной задачи.

Относительно системы уравнений отметим два обстоятельства.

  • 1. Система имеет очень высокую размерность (М - 1) х (N - 1), а традиционные методы точного решения (например, метод Гаусса) требуют для решения число алгебраических операций, пропорциональное третьей степени размерности системы.
  • 2. В матрице системы много нулевых элементов (неплотная матрица). Это обстоятельства позволяет разработать экономичные методы приближенного решения.

Рассмотренная постановка разностной задачи характерна для эллиптических уравнений. В газовой динамике такой вид имеет уравнение для функции тока или для потенциала скорости. В других разделах мы рассмотрим эффективные методы разрешения таких разностных схем.


Рис. 2.8.

П р и м с р 2. Разностная схема для простейшего параболического уравнения (нестационарная теплопроводность в стержне единичной длины).

Рассмотрим следующую задачу:


Отмстим, что в случае параболического уравнения имеем открытую область. При построении разностной схемы возникают несколько вариантов связи между разностными производными по пространству и по времени.

Проинтегрируем уравнение в пределах одного временного шага:


В зависимости от того, какую квадратурную формулу мы примем для вычисления интеграла в правой части, получим разные разностные схемы (рис. 2.9).

Связывая разностную производную по времени с пространственной производной, определенной на п -м временном слое, получим

явную ’разностную схему

Это эквивалентно приближенному вычислению интеграла, стоящего в правой части (2.12), но способу левых прямоугольников.


Рис. 2.9. Сетка и шаблоны для уравнения теплопроводности: а - область и сетка; б - шаблон явной схемы; в - шаблон неявной схемы; г - шаблон семейства шеститочечных схем; д - шаблон схемы

«чехарда»

В приведенной формуле заключен и метод решения сеточных уравнений:

Значение сеточной функции на следующем временном слое

определяется через известные значения гф на предыдущем. Перемещаясь последовательно слоями от начального условия и(х , 0) = у(х), можно найти решение во всей расчетной области. Разностный шаблон для этой схемы приведен на рис. 2.9, б .

Оценивая интеграл через значение подынтегральной функции па слое п + 1, используем разностный шаблон вида рис. 2.9, б, а разностный аналог дифференциального уравнения примет вид

Для того чтобы найти значения сеточной функции на следующем временном слое, при использовании этой разностной схемы необходимо решить совместно столько уравнений вида (2.14), сколько внутренних узлов расположено на п - 1-1 -м временном слое. С учетом краевых условий = / п+1 , Мд Г +1 = m n+1 система позволяет построить решение на следующем временном слое при известных значениях сеточной функции на предыдущем. Передвигаясь от начальных значений слоями, на каждом из которых необходимо решать систему уравнений, можно построить приближенное решение во всей области.

Рассмотренная разностная схема представляет собой пример неявной разностной схемы, ее называют схемой с опережением или чисто неявной схемой.

Шеститочечный разностный шаблон порождает семейство разностных схем, частными случаями которого являются две предыдущие:


При а = 0 имеем явную схему, при а = I - неявную с опережением, при а > 0 - неявную. При а - 0,5 получаем широко известную в вычислительной практике симметричную схему Кранка Николсона.

Приведенные схемы, разумеется, не исчерпывают всего многообразия разностных схем, основанных на разностной аппроксимации дифференциальных операторов. Вот пример явной разностной схемы, основанной на центрировании временной производной, схемы, использующей сеточную функцию на трех временных слоях:

Разностный шаблон захватывает три временных слоя. Схема имеет второй порядок аппроксимации как по времени, так и по пространственной переменной и является явной. Эта схема обладает рядом существенных недостатков, от большей части которых можно избавиться, заменив и ” в аппроксимации пространственной производной средним значением по двум временным слоям:

Полученная таким образом явная трехслойная схема

называется схемой Дюфортпа - Франкела , а отсутствие в ней значения сеточной функции в центральном узле объясняет название «чехарда», которое иногда применяется для схем такого рода.

На примерах было показано, что для одной и той же краевой задачи можно записать несколько разных разностных схем, т.е. в распоряжении исследователя имеется достаточно большой их выбор. Каким же условиям должна удовлетворять разностная схема, чтобы разностное решение соответствовало решению исходной дифференциальной задачи? Этот вопрос будет рассмотрен в следующем разделе.

Раздел ¹ 10. Численное решение уравнений в частных производных

Разностные схемы для уравнений эллиптического типа

Различные краевые задачи и аппроксимация граничных условий

Построение разностной схемы в случае задачи Дирихле для уравнения Пуассона

Метод матричной прогонки

Итерационный метод решения разностной схемы для задачи Дирихле

Уравнение параболического типа. Явные и неявные конечноразностные методы

Методы прогонки для уравнения параболического типа

Предметный указатель

Разностные схемы. Основные понятия

Пусть Д - некоторая область изменения независимых переменных x, y, ограниченная контуром. Говорят, что в области Д задано линейное дифференциальное уравнение второго порядка для функции U(x, y), если для любой точки из области Д имеет место соотношение

∂2 U

∂2 U

∂2 U

∂x2

∂x2

G(x, y)U = f(x, y),

где a(x, y), b(x, y), . . . - коэффициенты, f(x, y) - свободный член уравнения. Эти функции известны и их обычно считают определенными в замкнутой области Д = Д + .

График решения представляет собой поверхность в пространстве Oxyz.

Назад Первая Предыдущая Следующая Последняя Перейти Предметный указатель

Обозначим δ(x, y) = b2 − ac. Уравнение L(U) = f называется эллиптическим, параболическим или

гиперболическим в Д, если соответственно выполняются условия δ(x, y) < 0, δ(x, y) = 0, δ(x, y) > 0 для

всех (x, y) Д.

В зависимости от типа дифференциального уравнения по-разному ставятся граничные начальные

(10.1):

Уравнение Пуассона (уравнение эллиптического типа)

∂2 U ∂2 U

∂x 2 + ∂y 2 = f(x, y)

Назад Первая Предыдущая Следующая Последняя Перейти Предметный указатель

Уравнение теплопроводности (уравнение параболическго типа)

∂U = ∂ 2 U + f(x, t) ∂t ∂x2

Волновое уравнение (уравнение гиперболического типа)

∂2 U ∂2 U

∂x 2 + ∂y 2 = f(x, y)

Сходимость, аппроксимация и устойчивость разностных схем

Пусть U есть решение дифференциального уравнения

заданного в Д. Рассмотрим некоторое множество Дh = {Mh } состоящее из изолированных точек Mh , принадлежащих замкнутой области Д = Д + . Число точек в Дh , будем характеризовать величиной h; чем меньше h, тем большим будет число точек в Дh . Множество Дh называется сеткой, а точки Mh Дh - узлами сетки. Функция, определенная в узлах, называется сеточной функцией. Обозначим через U пространство непрерывных в D функций V (x, y). Через Uh обозначим пространство, образованное совокупностью сеточных функций Vh (x, y), определенных на Дh . В методе сеток осуществляется замена пространства U на пространство Uh .

Пусть U(x, y) - точное решение уравнения ((10.2 )) и U(x, y) принадлежит U. Поставим задачу отыскания значений Uh (x, y). Эти значения в совокупности образуют таблицу, в которой число значений

Назад Первая Предыдущая Следующая Последняя Перейти Предметный указатель

равно числу точек в Дh . Точно поставленную задачу удается решить редко. Как правило, можно вычислить некоторые сеточные значения U(h) , относительно которых можно предполагать, что

U(h) ≈ Uh (x, y).

Величины U(h) называются приближенными сеточными значениями решения U(x, y). Для их вычисления строят систему численных уравнений, которую мы будем записывать в виде

Lh (U(h) ) = fh ,

есть разностный оператор,

соответствующий оператору

зуется по F аналогично тому, как U

образовывалось по U. Формулу (10.3 ) будем называть разностной

схемой. Пусть в линейных пространствах Uh и Fh введены соответственно нормы k · kU h и k · kF h , которые являются сеточными аналогами норм k · kU и k · kF в исходных пространствах. Будем говорить, что разностная схема (10.3 ) является сходящейся, если при h → 0 выполняется условие

kUh (x, y) − Uh kU h → 0.

Если выполняется условие

kUh (x, y) − Uh kU h 6 chs ,

где c - постоянная, не зависящая от h и s > 0, то говорят, что имеет место сходимость со скоростью порядка s относительно h.

Говорят, что разностная схема (10.3 ) аппроксимирует задачу (10.2 ) на решении U(x, y), если

Lh (Uh (x, y)) = f(h) + δf(h) и

δf(h) F h → 0 приh → 0.

Назад Первая Предыдущая Следующая Последняя Перейти Предметный указатель

Величина δf(h) называется погрешностью аппроксимации или невязкой разностной схемы. Если

δf (h) F h 6 Mh σ, где M - константа, не зависящая от h и σ > 0, то говорят, что задана разностная схема (10.3 ) на решении U(x, y) с погрешностью порядка σ относительно h.

Разностная схема (3) называется устойчивой, если существует такое h0 > 0, что для всех h < h0 и любых f(h) Fh выполняются условия

Разностная схема (10.3 ) имеет единственное решение;

U (h) U h

f(h) F h , где M - постоянная, не зависящая от h и f(h) .

Иначе говоря, разностная схема является устойчивой, если ее решение непрерывно зависит от входных данных. Устойчивость характеризует чувствительность схемы к различного рода погрешностям, она является внутренним свойством разностной задачи и это свойство не связывается непосредственно с исходной дифференциальной задачей, в отличие от сходимости и аппроксимации. Между понятиями сходимости, аппроксимации и устойчивости существует связь. Она состоит в том, что из аппроксимации и устойчивости следует сходимость.

Теорема 1 Пусть разностная схема L h (U h (x, y)) = f (h) аппроксимирует задачу L(U) = f на решении U(x, y) с порядком s относительно h и устойчива. Тогда эта схема будет сходиться, и порядок ее сходимости будет совпадать с порядком аппроксимации, т.е. будет справедлива оценка

Uh (x, y) − Uh U h 6 khs ,

где k - постоянная, не зависящая от h .

Доказательство . По определению аппроксимации имеем

(h) F h 6 M(Chs ) = Khs ,

где K = MC. Таким образом, оценка (10.4 ) установлена и теорема доказана. Обычно применение метода сеток заключается в следующем:

1. Вначале указывается правило выбора сетки, т.е. указывается метод замены области Д и контура Г некоторой сеточной областью. Чаще всего сетка выбирается прямоугольной и равномерной.

2. Затем указывается и строится конкретно одна или несколько разностных схем. Проверяется условие аппроксимации и устанавливается ее порядок.

3. Доказывается устойчивость построенных разностных схем. Это один из наиболее важных и сложных вопросов. Если разностная схема обладает аппроксимацией и устойчивостью, то о сходимости судят по доказанной теореме.

4. Рассматривается вопрос численного решения разностных схем.

В случае линейных разностных схем это будет система линейных алгебраических уравнений. Порядок таких систем может быть большим.

Назад Первая Предыдущая Следующая Последняя Перейти Предметный указатель

Сетка и шаблон. Для большинства разностных схем узлы сетки лежат на пересечении некоторых прямых линий (в многомерных задачах – гиперплоскостей), проведенных либо в естественной системе координат, либо в специально подобранной по форме области G .

Если одна из переменных имеет физический смысл времени t , то сетку обычно строят так, чтобы среди ее линий (или гиперплоскостей) были линии t = t m . Совокупность узлов сетки, лежащих на такой линии или гиперплоскости, называют слоем.

На каждом слое выделяют направления, вдоль которых меняется только одна пространственная координата. Например, для переменных x , y , t есть направления x (t = const , y = const ) и направление y (t = const , х = const ).

Составляя разностные схемы (26.2) и (26.4), мы использовали во всех внутренних узлах области однотипную разностную аппроксимацию производных. Иными словами, при написании каждого разностного уравнения около некоторого узла сетки бралось одно и то же количество узлов, образующее строго определенную конфигурацию, которую мы назвали шаблоном данной разностной схемы (см. рис. 26.2).

Определение. Узлы, в которых разностная схема записана на шаблоне, называются регулярными, а остальные – нерегулярными.

Нерегулярными являются обычно граничные узлы, а иногда также лежащие вблизи границы узлы (такие, что взятый около этого узла шаблон выходит за границу области).

Составление разностной схемы начинается с выбора шаблона. Шаблон не всегда определяет разностную схему однозначно, но существенно влияет на ее свойства; например, далее мы увидим, что на шаблоне рис. 26.2b нельзя составить хорошей разностной схемы для задачи теплопроводности (26.1). Для каждого типа уравнений и краевых задач требуется свой шаблон.

Явные и неявные разностные схемы

Обсудим вопрос о фактическом вычислении разностного решения. Большая часть физических проблем приводит к уравнениям, содержащим время в качестве одной из переменных. Для таких уравнений ставится обычно смешанная краевая задача, типичным случаем которой является задача теплопроводности (26.1).

К подобным задачам применяют послойный алгоритм вычислений. Рассмотрим его на примере схем (26.2) и (26.4).

В схеме (26.4) на исходном слоеm = 0 решение известно в силу начального условия. Положим m = 0 в уравнениях (26.4). Тогда при каждом значении индекса n уравнение содержит одно неизвестное ; отсюда можно определитьпри
Значенияиопределяются по краевым условиям (26.3). Таким образом, значения на первом слое вычислены. По ним аналогичным образом вычисляется решение на втором слое и т.д.

Схема (26.4) в каждом уравнении содержит только одно значение функции на следующем слое; это значение нетрудно явно выразить через известные значения функции на исходном слое, поэтому такие схемы называются явными.

Схема (26.2) содержит в каждом уравнении несколько неизвестных значений функции на новом слое; подобные схемы называются неявными. Для фактического вычисления решения перепишем схему (26.2) с учетом краевого условия (26.3) в следующей форме

(26.5)

На каждом слое схема (26.5) представляет собой систему линейных уравнений для определения величин
; правые части этих уравнений известны, поскольку содержат значения решения с предыдущего слоя. Матрица линейной системы трехдиагональная, и решение можно вычислить алгебраической прогонкой.

Рассмотренный сейчас алгоритм достаточно типичен. Он используется во многих неявных разностных схемах для одномерных и многомерных задач. Дальше мы будем вместо индекса m часто применять сокращенные обозначения

В этих обозначениях явная и неявная разностные схемы принимают соответственно следующий вид


Невязка. Рассмотрим операторное дифференциальное уравнение общего вида (не обязательно линейное)

Au = f , или Au f = 0.

Заменяя оператор А разностным оператором A h , правую часть f – некоторой сеточной функцией , а точное решениеu – разностным решением y , запишем разностную схему

или
. (26.6)

Если подставить точное решение u в соотношение (26.6), то решение, вообще говоря, не будет удовлетворять этому соотношению
. Величину

называют невязкой.

Невязку обычно оценивают при помощи разложения в ряд Тейлора. Например, найдем невязку явной разностной схемы (26.4) для уравнения теплопроводности (26.1а). Запишем это уравнение в каноническом виде

Поскольку в данном случае
то

Разложим решение по формуле Тейлора около узла (x n , t m ), предполагая существование непрерывных четвертых производных по х и вторых по t

(26.7)

где

Подставляя эти разложения в выражение невязки и пренебрегая, в силу непрерывности производных, отличием величин
от (x n , t m ) найдем

(26.8)

Таким образом, невязка (26.8) стремится к нулю при
и
Близость разностной схемы к исходной задаче определяется по величине невязки. Если невязка стремится к нулю приh и стремящихся к нулю, то говорят что такая разностная схема аппроксимирует дифференциальную задачу. Аппроксимация имеетр -й порядок, если
.

Выражение (26.8) дает невязку только в регулярных узлах сетки. Сравнивая (26.3) и (26.1б), легко найдем невязку в нерегулярных узлах

Замечание 1. Решение задачи теплопроводности с постоянным коэффициентом (26.1) в области непрерывно дифференцируемо бесконечное число раз. Однако учет пятых и более производных в разложении в ряд Тейлора (26.7) прибавит к невязке (26.8) только члены более высокого порядка малости поиh , т.е. по существу, не изменит вида невязки.

Замечание 2. Пусть по каким-либо причинам решение исходной задачи дифференцируемо небольшое число раз; например, в задачах с переменным коэффициентом теплопроводности, гладким, но не имеющим второй производной, решение имеет лишь третьи непрерывные производные. Тогда в разложении в ряд Тейлора (26.7) последними будут члены
не точно компенсирующие друг друга. Это приведет к появлению в невязке (26.8) члена типа
т.е. невязка будет иметь меньший порядок малости, чем для четырежды непрерывно дифференцируемых решений.

Замечание 3. Преобразовав выражение невязки с учетом того, что входящая в него функция u (x ,t ) есть точное решение исходного уравнения и для нее выполняются соотношения

Подставляя это выражение в (26.8), получим

Если выбрать шаги по пространству и времени так, чтобы
то главный член невязки обратится в нуль и останутся только члены более высокого порядка малости поиh (которые мы опускали). Этим приемом пользуются при построении разностных схем повышенной точности.