Методы составления разностных схем

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

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

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

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

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

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

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

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

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

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

. (35)

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

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

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

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

(36)

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

(37)

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

, (38)

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

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

Листинг_№4

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

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

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

function ktmp

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);

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

Tm=7;

for i=1:N

if x(i)<=0.5*a

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

end

if x(i)>0.5*a

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

end

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

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

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

hold on

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

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

for t=1:20

for n=2:(N-1)

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);

end

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

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

for n=2:(N-1)

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));

end

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

y(N)=0;

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

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

end

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

plot(x,y);

hold on

end

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

function z=k(x)

global a k1 k2 k3

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

z=k1;

end

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

z=k2;

end

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

z=k3;

end

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

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

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

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

. (39)

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

. (40)

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

(41)

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

т.е. при

. (42)

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

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

. (43)

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

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

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

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

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

. (44)

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

. (45)

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

(46)

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

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

(47)

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

. (48)

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

(49)

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

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

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

(50)

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

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