Явная разностная схема для уравнения теплопроводности. Разностная схема Метод золотого сечения

Используя шаблон для каждого внутреннего узла области решения апроксимируется уравнение теплопроводности

Отсюда найдем:

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

Затем с помощью соотношений

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

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

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

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

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

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

Рассмотрим разностную схему, использующую другой шаблон

Метод 36

Неявная разностная схема для уравнения теплопроводности.

Подставим в уравнение теплопроводности:

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

Схема решения задачи следующая:

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

Отличие от явной схемы - значения на очередном временном уровне вычисляются не непосредственно с помощью готовой формулы, а находится путем решения системы уравнений, т.е. значения неизвестных находятся неявно путем решения СЛАУ. Поэтому разностная схема называется неявной. В отличие от явной неявная является абсолютно устойчивой.

Тема №9

Задачи оптимизации.

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

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

- целевая функция.

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

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

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

Теория и методы решения задач оптимизации при наличии ограничений составляют предмет исследований одного из разделов прикладной математики – математического программирования.

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

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

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

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

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

То есть, делаем пока выполняется (*)

При этом производится вычислений целевой функции.

Требуется найти такое значение, чтобы выражение (*) было получено при наименьшем

числе вычислений .

Метод 37

Метод половинного деления.

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

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

Метод 38

Метод золотого сечения.

Одним из эффективных способов является метод золотого сечения. Золотым сечением отрезка называется точка для которой выполняется условие

Таких точек две: =0,382 +0,618

0,618 +0,382 .

Отрезок делится точками и а после находится точка, целевая функция в которой максимальна. В результате чего находится изменённый отрезок длинною 0,618( - ) .

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

Метод 39

Метод покоординатного подъёма (спуска).

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

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

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

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

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

Методом разностной аппроксимации мы уже пользовались при составлении схем (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 рассмотрим решение , соответствующее возмущенной правой части . В дальнейшем будем считать, что . Это можно предположить, т.к. исследуется устойчивость по правой части.

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

ГЛАВА 4. ЭЛЕМЕНТАРНЫЕ ПРИМЕРЫ РАЗНОСТНЫХ СХЕМ

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

§ 8. Понятие о порядке точности и об аппроксимации

1. Порядок точности разностной схемы.

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

Начнем с простейшей разностной схемы, основанной на использовании разностного уравнения

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

откуда находим решение уравнения (2) при начальном условии :

Точное же решение задачи (1) имеет вид . Оно принимает в точке значение

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

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

Таким образом, равенство (3) примет вид

т. е. погрешность (5) стремится к нулю при и величина погрешности имеет порядок первой степени шага.

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

Решим теперь задачу (1) с помощью разностного уравнения

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

Не ясно, как задавать их. Чтобы разобраться в этом, воспользуемся явной формой решения уравнения (7) (см. § 3 формулы ):

Разложения (9) по формуле Тейлора корней характеристического уравнения позволяют дать приближенные представления для Проведем подробно вывод такого представления -

Так как , то

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

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

Все дальнейшие выводы мы будем получать путем исследования этой формулы.

Заметим, что если коэффициент стремится при к конечному пределу b, то первое слагаемое правой части равенства (12) стремится к искомому решению задачи (1).

Сетка и шаблон. Для большинства разностных схем узлы сетки лежат на пересечении некоторых прямых линий (в многомерных задачах – гиперплоскостей), проведенных либо в естественной системе координат, либо в специально подобранной по форме области 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 (которые мы опускали). Этим приемом пользуются при построении разностных схем повышенной точности.

1. В системе координат xOt строим прямоугольную сетку с шагом h по оси Ох и с шагом τ по оси Ot :

a) x i =ih , i = l, n , n=L/h ;

б) t k =k τ, k= l, m , m =T/τ;

в) и i , k = u (x i , t k ) = u (ih , k τ).

2. Вычисляем значения функции u (x i , t k ) в узлах, лежащих на прямых х= 0 и x=L :

3. Вычисляем u i ,0 =f (ih ), i= 1, n .

4. Используя (1.16) или (1.23), найдем решение для всех внутренних узлов: u i , k + n , i= l, n -l , k= 0, m -l .

1.3. Решение смешанной задачи для волнового уравнения методом сеток

1.3.1. Постановка задачи. Алгоритм метода

Рассмотрим смешанную задачу (т. е. заданы начальные и граничные условия) для волнового уравнения

в области D ={0≤x≤L , 0≤t≤T } с начальными условиями

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

Будем предполагать, что f (x ), g (x ) – достаточно гладкие функции, причем выполнены условия согласования в двух углах области D (x =0, t =0), (x=L , t =0), обеспечивающие существованне и единственность решения u (x , t ).

Для дискретизации исходной задачи построим в области

прямоугольную сетку

где h – шаг сетки в направлении х , τ – шаг сетки в направлении t ,

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

которые аппроксимируют волновое уравнение (1.24) в узле (х i , t k ) с погрешностью O (h 2 + τ 2).

Здесь u i , k – приближенное значение функции и (х , t ) в узле (x i , t k ).

Полагая λ = аτ/h , получим трехслойную разностнуюсхему:

Схема (1.28) называется трехслойной потому, что связывает между собой значения u i , k функции и (х , t ) на трех временных слоях с номерами (k -l), k , (k +1).

Разностной схеме (1.28) соответствует пятиточечный трехслойный шаблон типа «крест» (рис. 1.2).

Схема (1.28) связывает значения u i , k =u (ih , kτ ) на трех слоях по времени, и чтобы перейти на уровень (k +1), необходимо знать как u i , k , так и u i , k -1 , что является следствием того, что дифференциальное уравнение (1.24) содержит вторую производную по времени. Численное решение задачи (1.24) – (1.26) состоит в вычислении приближенных значений u i , k решения u (х , t ) в узлах (х i , t ) при i = 1, n , k =1, m . Схема вычислений по (1.28) является явной, она позволяет вычислить приближенно значения функции в узлах (k +1)-го слоя по известным ее значениям на двух предыдущих слоях. На первых двух слоях значения функции определяются из начальных условий (1.25). Полагаем

Для производной по времени применяем аппроксимацию (1.5)

Порядок аппроксимации (1.30) равен О (τ).

Заметим, что (1.29), (1.31) дают решения для первыхдвухстрок: k =0, k =1. Подставляя k= 1 в (1.28), получим:

Все слагаемые в правой части уравнения (1.32) включают значения и i , k только из первых двух строк сетки; но ведь все эти значения известны из начальных условий.

После этого, зная решения и i ,1 , и i ,2 , можно по (1.28) вычислить значения функции и i , k на третьем временном слое, четвертом и т. д.

Описанная выше схема вычислений (1.28) – (1.31) aппpoксимирует задачу (1.24) – (1.26) с точностью О (τ+h 2). Невысокий порядок аппроксимации по τ объясняется использованием слишком грубой аппроксимации для производной по t в формуле (1. 30).

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

Это означает, что при выполнении (1.33) малые погрешности, возникающие, например, при вычислении на первом слое, не будут неограниченно возрастать при переходе к каждому новому временному слою. При выполнении условия Куранта разностная схема (1.28) обладает равномерной сходимостью, т. е. при h →0 и τ→0 решение разностной задачи (1.28) – (1.31) равномерно стремится к решению исходной задачи (1.24) – (1.26).

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

Таким образом, как только выбрана величина шага h в направлении х , то появляется ограничение на величину шага τ по времени. Отличительная особенность всех явных методов заключается в том, что при их использовании должно соблюдаться некоторое условие типа (1.33), обеспечивающее сходимость и устойчивость метода.