Конечные разности и разностные уравнения

КОНЕЧНЫЕ РАЗНОСТИ

КОНЕЧНЫЕ РАЗНОСТИ. Исчисление конечных разностей связано с изучением свойств и применений разностей между соседними членами какой-нибудь последовательности или между значениями функции в точках, расположенных с постоянным интервалом в некотором пространстве. Слово «конечные» используется здесь в несколько устаревшем смысле «не бесконечно малые», т.е. не связанные с предельными переходами. Поскольку дифференциальное исчисление занимается изучением пределов разностей, а исчисление конечных разностей – самими разностями, то естественно, что между этими двумя теориями существуют много параллелей. Исчисления конечных разностей используются при интерполяции в математических таблицах, при суммировании числовых рядов, при вычислении интегралов и дифференцировании функций. Разности встречаются также в любой ситуации, когда надо описать поведение объекта, который испытывает воздействие меняющихся условий на определенном расстоянии (во времени и в пространстве). Например, термостату требуется значительное время, чтобы отреагировать на изменение температуры, поэтому он реагирует не на текущую температуру, а на ту, что была минуту назад. Другой пример: автомашиной управляет водитель, которому требуется какое-то время, чтобы отреагировать на возникшую на дороге ситуацию.

Под конечной разностью первого порядка функции f (x) принято понимать величину

где d – некоторая постоянная, которую часто, но не всегда, принимают равной 1. Разность второго порядка обозначается D2f и представляет собой разность разностей, т.е.

Продолжив этот процесс, мы получим разности более высоких порядков D3f (x), D4f (x), ј .

Данные выше определения можно также применить к членам любых последовательностей величин, например, к последовательности

3, 6, 11, 18, 27, 38, ј .

Первые разности равны

6 – 3, 11 – 6, 18 – 11, 27 – 18, 38 – 27, ј,

разности второго порядка постоянны и равны 2.

В общем виде такие последовательности можно записать как

где разности первого, второго и т.д. порядков определяются выражениями

а n может принимать любое допустимое для индекса значение.

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

где индексы могут принимать любые убывающие значения. В этом случае вместо символа D используется символ «разделенной разности». Разделенные разности первого и второго порядков определяются следующим образом:

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

У истоков теории.

Хотя исследование свойств и использование конечных разностей приходится на современный период развития математики, Птолемей (ок. 150 н.э.) ввел в Альмагесте таблицу разностей первого порядка, чтобы облегчить расчеты в таблице длин хорд. Разности второго порядка использовал при вычислении своих таблиц логарифмов в 1624 Г.Бриггс. Теория интерполяции берет начало со знаменитой пятой леммы из 3-й книги Математических начал (1687) И.Ньютона, в которой впервые была приведена формула, носящая ныне его имя. Частный случай формулы Ньютона, открытый также независимо его современником Дж.Грегори (1638–1675), приведен ниже (см. формулу (7)). В общей формуле интерполяции Ньютона использовались разделенные разности, хотя этот термин, по-видимому, был введен О.де Морганом (1806–1871) в 1848. Первое применение исчисления конечных разностей к задачам теории вероятностей принято связывать с именами П.де Монтмора (1678–1719) и А.де Муавра (1667–1754).

Хотя Л.Эйлер (1707–1783) в своих работах по дифференциальному исчислению использовал предельные переходы в конечных разностях, основания современной теории конечных разностей были заложены в основном Ж.Лагранжем (1736–1813) и П.Лапласом (1749–1827). Первый из них ввел в исчисление конечных разностей символические методы, второй сделал конечные разности главным инструментом в своей Аналитической теории вероятностей (1812).

Под влиянием этих работ математики 19 в. принялись интенсивно разрабатывать предмет, и в 1860 Дж.Буль выпустил свой классический Трактат об исчислении конечных разностей. С тех пор это исчисление и круг его приложений существенно расширились. Одно из наиболее важных приложений конечные разности нашли в статистике. Особенно полезными они оказались в теории сериальной корреляции, в анализе случайных последовательностей и статистических временных рядов.

Интерполяция.

Чтобы понять, как конечные разности используются при интерполяции, рассмотрим следующую таблицу:

Методы конечных разностей

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

Рассмотрим сущность такого метода решения для дифференциального уравнения второго порядка (1.41) при заданных граничных условиях (1.42). Разобьем отрезок [0,1] на n равных частей точками xi= ih(i= 0,1. , n). Решение краевой задачи (1.41), (1.42) сведем к вычислению значений сеточной функции yi в узловых точках xi. Для этого напишем уравнение (1.42) для внутренних узлов:

(1.49)

Заменим производные, входящие в эти соотношения, их конечно-разностными аппроксимациями:

(1.50)

Подставляя эти выражения в (1.49), получаем систему разностных уравнений:

(1.51)

являющуюся системой n-1 алгебраических уравнений относительно значений сеточной функции Входящие в данную систему y0 (при i = 1) и уп (при i = п — 1) берут из граничных условий (1.42):

На практике часто граничные условия задают в более общем виде (1.38):

(1.52)

В этом случае граничные условия также должны представляться в разностном виде путем аппроксимации производных Y‘(0) и Y‘(1) с помощью конечно-разностных соотношений. Если использовать односторонние разности (соответствующий шаблон показан на рис. 1.7, а),при которых производные аппроксимируются с первым порядком точности, то разностные граничные условия примут вид

(1.53)

Из этих соотношений легко находятся значения y0, yn.

Однако, как правило, предпочтительнее аппроксимировать производные, входящие в (1.52), со вторым порядком точности с помощью центральных разностей:

Рис. 1.7. Аппроксимация граничных условий

В данные выражения входят значения сеточной функции и yn+1 в так называемых фиктивных узлах х=1-hи х =1+h, лежащих вне рассматриваемого отрезка (рис. 1.7, б). В этих узлах значения искомой функции также должны быть найдены. Следовательно, количество неизвестных значений сеточной функции увеличивается на два. Для замыкания системы привлекают еще два разностных уравнения (1.51) при i = 0, i = п.

Аппроксимировать граничные условия со вторым порядком можно и иначе (см. рис. 1.7, в). В этом случае используют аппроксимации:

Таким образом, решение краевой задачи для дифференциального уравнения сведено к решению системы алгебраических уравнений вида (1.51). Эта система является линейной или нелинейной в зависимости от того, линейно или нелинейно дифференциальное уравнение (1.41). Методы решения таких систем рассмотрены ранее.

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

(1.54)

с граничными условиями вида

Разобьем отрезок [0,1] на части с постоянным шагом h с помощью узлов . Аппроксимируем вторую производную Y²конечно-разностным соотношением (1.50). При этом значения искомой функции в узлах Y(xi) приближенно заменяем соответствующими значениями сеточной функции yi. Записывая уравнение (1.54) в каждом узле с использованием указанных аппроксимаций, получаем

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

(1.56)

Получилась система n—1 линейных уравнений, число которых совпадает с числом неизвестных значений сеточной функции в узлах. Ее значения на концах отрезка определены граничными условиями (1.55):

Решив систему уравнений (1.56) с учетом условий (1.57), найдем значения сеточной функции, которые приближенно равны значениям искомой функции. Покажем, что такое решение существует и сходится к точному решению при h→ 0.

Для доказательства существования решения рассмотрим систему линейных уравнений (1.56). Ее матрица является трехдиагональной; на главной диагонали находятся элементы . Поскольку р(х) > 0, то pi > 0, и диагональные элементы матрицы преобладают над остальными, так как в каждой строке модули этих элементов больше суммы модулей двух остальных элементов, каждый из которых равен единице. При выполнении этого условия решение системы линейных уравнений существует и единственно.

Что касается сходимости решения, то здесь имеет место следующее утверждение.

Утверждение. Если функции р(х) и f(x) дважды непрерывно дифференцируемы, то при h→0 разностное решение равномерно сходится к точному со скоростью O(h2).

Это — достаточное условие сходимости метода конечных разностей для краевой задачи (1.54), (1.55).

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

Этот метод на практике используется также и при р(х) Будет полезно почитать по теме:

Конечно-разностные аппроксимации производных (стр. 1 )

Из за большого объема этот материал размещен на нескольких страницах: 1 2 3 4 5 6 7

Конечно-разностные аппроксимации производных

Конечно-разностные аппроксимации производных (конечные разности) — способ приближенного вычисления частных производных

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

Или более коротко с использованием индексов точек:

Отсюда , где — остаток.

Отбрасывая остаток можно получить правую разность:

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

Аналогичным образом, разлагая в ряд функцию можно получить:

Получим новую аппроксимацию первой производной:

которая называется левой разностью. У нее погрешность также определяется членом, содержащим в первой степени. Однако, если из выражения (1) вычесть (2), то можно получить более точную аппроксимацию первой производной, которая называется центральной разностью:

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

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

Отброшенный остаток будет содержать член с во второй степени (после деления на )

Исходя из определения, можно получить выражения для третьей, четвертой и более старших разностей:

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

правая разность по оси : ; правая разность по оси : ; левая разность по оси : ; левая разность по оси : ; центральная разность по оси : ; центральная разность по оси : ; вторая разность по оси : ; вторая разность по оси : .

Смешанная производная может быть получена следующим образом:

Алгоритм решения стационарных краевых задач методом конечных разностей

Метод конечных разностей — универсальный сеточный численный метод решения задач микроуровня.

Алгоритм решения стационарных краевых задач методом конечных разностей — последовательность действий, приводящая к решению стационарной задачи микроуровня

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

Решение линейных одномерных стационарных краевых задач с помощью МКР

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

Одномерное стационарное уравнение теплопроводности для изотропной среды выглядит следующим образом:

В соотвествии с алгоритмом решения стационарных краевых задач методом конечных разностей наносим на объект равномерную сетку, как это показано на рис. 1.

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

для узла 1:
для узла 2:

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

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

Рассмотрим теперь решение задачи с краевым условием второго рода, на правой границе стержня задан тепловой поток:

Пусть и .

Запишем разностные аналоги для внутренних узлов сетки:

для узла 1:
для узла 2

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

Проще всего воспользоваться левой разностью:

Решая эту систему уравнений, получим , , .

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

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

Теперь можно записать следующую систему конечно-разностных уравнений:

для узла 1:
для узла 2:
для узла 3:
граничное условие второго рода:

За повышение точности пришлось заплатить увеличением размерности системы конечно-разностных уравнений.

Решение нелинейных одномерных стационарных краевых задач с помощью МКР

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

Одномерное стационарное уравнение теплопроводности для анизотропной среды выглядит следующим образом:

где — коэффициент теплопроводности.

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

Предположим, что толщина стекла . Температура в комнате , на улице —

Тепловой поток на улицу пропорционален градиенту температуры, то есть .

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

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

для узла 1:
для узла 2:
для узла 3:

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

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

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

для узла 1:
для узла 2:

Данную систему придется решать итерационными методами.

Решение нестационарных одномерных задач с помощью МКР

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

Но поскольку при этом возможны различные сочетания конечных разностей по оси координат и по времени, можно получить различные схемы решения нестационарных задач. Рассмотрим их на примере нестационарного уравнения теплопроводности:

Пусть при записи разностей нижний индекс соответствует оси , а верхний — оси времени.

Первый вариант разностного уравнения, апроксимирующего исходное (1):

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

Второй вариант разностного уравнения, апроксимирующего исходное (1):

называется неявной разностной схемой, поскольку в этом уравнении несколько неизвестных величин, относящихся к -му временному слою. Для их нахождения придется записать систему разностных уравнений для всех внутренних узлов сетки, и решить ее.

Графическое изображение разностных уравнений получило название шаблонов решения сответствующих задач. В данном случае на рис. 1,а представлен шаблон явной разностной схемы, а на рис. 1,б — неявной.

Рис. 1. Шаблоны явной и неявной разностной схемы

Использование шеститочечного шаблона применено в схеме Кранка-Николсона:

В общем случае использования шеститочечного шаблона, имеем схему с весами:

которая при является неявной.

Примеры решения нестационарных задач с помощью МКР

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

Решим задачу с помощью явной разностной схемы.

Одномерное нестационарное уравнение теплопроводности для изотропной среды выглядит следующим образом:

Пусть , выберем значения шага по оси и значение шага по оси времени .

Наносим на объект равномерную сетку по оси , как это показано на рис. 1.

Записываем явную разностную схему для узла 1:

где — граничное условие, — начальные условия, отсюда .

Записываем явную разностную схему для узла 2:

где — граничное условие, — начальные условия, отсюда .

Таким образом найдено температурное поле в момент времени .

Аналогично для момента времени :

Для момента времени :

Получили картину прогревания стержня в течение трех единиц времени, представленную на рис. 2.

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

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

для момента времени ; для момента времени для момента времени .

Теперь картина прогревания не противоречит физическому смыслу задачи.

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

вычислить значения производных по времени во всех внутренних узлах объекта; определить максимальное из этих значений; разрешить измениться переменной в этом узле на некоторую заданную величину, которая определяется из физического смысла задачи. (Например для нашей задачи максимальной значение температуры внутри стержня , за один шаг по времени можем позволить измениться ей, допустим, на . Исходя из этого вычисляем значение ); выполняем шаг по времени для всех узлов, изменение температуры во всех узлах не превысит разрешенной величины; если модельное время не закончилось переходим к пункту 1.

Рассмотрим решение задачи явной разностной схемой с граничными условиями второго рода (типа Неймана).

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