Прямые методы решения

Стационарные (не эволюционные) схемы можно составлять непосредственно, аппроксимируя производные разностями, а также с помощью интегро-интерполяционного метода. Так для многомерного уравнения D u = - f простейшая разностная схема имеет следующий вид:

. (43)

Для стационарных схем типа (43) наиболее сложным является вопрос о реальном вычислении решения. Действительно, разностная схема (43) имеет N^p неизвестных, если считать, что число узлов по каждому из направлений N. На листинге_№6 приведен код программы рисования матрицы двумерного уравнения (43). В программе используется специальная функция spy пакета MATLAB, которая предназначена для изображения разреженных матриц.

Листинг_№6

%Программа рисования структуры матрицы

%двумерного разностного уравнения (43)

clear all

%Выбираем число узлов сеток по x1 и x2

N=7; M=7;

%Создаем матрицу двумерного разностного

%уравнения (43), заполненную нулями

A=zeros(N*M);

%Вносим единицы в те места матрицы, которые

%отвечают уравнениям разностной схемы (43)

%для регулярных узлов

for m=2:(M-1)

for n=2:(N-1)

k=n+N*(m-1);

kn1=n-1+N*(m-1);

kn2=n+1+N*(m-1);

km1=n+N*(m-2);

km2=n+N*m;

A(k,k)=1;

A(k,kn1)=1; A(k,kn2)=1;

A(k,km1)=1; A(k,km2)=1;

end

end

%Рисуем матрицу A в координатах n и m

spy(A);

%Создаем матрицу B заполненную нулями для

%учета нерегулярных узлов, т.е. для учета

%возможных граничных условий

B=zeros(N*M);

%Формируем матрицу B, добавляя в нее единицы,

%учитывающие возможные ненулевые граничные

%условия

for n=1:N

B(n,n)=1; B(n,n+N)=1;

B(n+N*(M-1),n+N*(M-1))=1;

B(n+N*(M-1),n+N*(M-2))=1;

end

for m=1:M

B(1+N*(m-1),1+N*(m-1))=1;

B(1+N*(m-1),2+N*(m-1))=1;

B(N+N*(m-1),N+N*(m-1))=1;

B(N+N*(m-1),N-1+N*(m-1))=1;

end

hold on

%Добавляем образ граничных условий B к

%матрице A

spy(B,10);

Результат работы кода программы листинга_№6 приведен на рис.6, где изображена матрица двумерного уравнения (43). Матрица ленточная, причем лента не полностью заполнена и имеет ширину ~ N. Красным цветом выделены элементы матрицы, наличие которых отвечает заданию возможных граничных условий второго рода при x ₁ = 0, x ₁ = a и x ₂ = 0, x ₂ = b.

Рис.6. Матрица двумерного разностного уравнения (43)

Вычисление разностного решения методом Гаусса, который не учитывает разреженность матрицы, требует ~ N ^{3 p} операций (при p = 1 на метод прогонки уходит ~ N операций), т.е. на каждый узел сетки приходится ~ N ^{2 p} операций. При заметном значении N это число операций может оказаться неприемлемо большим. Это означает, что прямое решение двумерного разностного уравнения (43) методом Гаусса возможно при значениях N, не более 100.

Если строить схемы высокого порядка точности, например, с помощью сплайнов в комбинации с вариационными методами и использовать сгущающиеся сетки (обычно в последовательности N = 4, 8, 16, 32) с уточнением решения по способу Рунге, то может быть достигнута заметная точность даже при небольшом числе узлов сетки.

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

Быстрое преобразование Фурье основано на том, что, если N_a — число узлов сетки по переменной x_a разлагается на множители, то вычислять коэффициенты дискретного преобразования Фурье можно не по прямым формулам, а по более экономичным, рекуррентным формулам. При метод является особенно быстрым и требует ~ 4log₂ N действий на каждый узел сетки.

Первоначально рассмотрим одномерную задачу для уравнения с постоянными коэффициентами . Аппроксимируем последнее уравнение на разностной схеме x_n = nh, n = 0,1,…, N:

. (44)

Разностное решение y_n будем искать в виде разложения Фурье:

, (45)

где . Подставим (45) в (44), умножим полученное уравнение на , просуммируем по n от 0 до N -1. Учитывая условие ортогональности гармоник, найдем

, (46)

где величины

(47)

являются дискретными коэффициентами Фурье правой части уравнения. Разностное решение уравнения (44) может быть получено согласно формулам (45) — (47). Формулы (45) — (47) являются неэкономичными, поскольку надо вычислить N коэффициентов b_p, на вычисление каждого из которых требуется ~ 2 N операций, т.е. всего ~ 2 N ² операций. Последний показатель заметно уступает методу прогонки, где требуется ~ N операций.

Пусть теперь число N можно разложить на сомножители, например, N = KL, тогда формулу (47) можно преобразовать таким образом, чтобы общее количество требуемых операций уменьшилось. Представим индексы n и p в виде:

(48)

С учетом представления (48) перепишем сумму (47) в виде двойной суммы:

. (49)

Поскольку , отбросим в (49) в показатели степени последнее слагаемое, тогда можно записать следующее представление

, (50)

где

. (51)

Проведем баланс числа операций по формулам (50), (51). Вычисление N коэффициентов b (p) по формуле (50) требует 2 NL операций. Вычисление KL = N вспомогательных коэффициентов по формуле (51) требует дополнительно 2 NK операций. Следовательно, всего операций для определения коэффициентов Фурье по формулам (50), (51) равно 2 N (L + K), что существенно меньше 2 N ². Действительно, пусть , тогда требуется операций, что в раз меньше, чем 2 N ².

Если K может быть также разложено на множители, то формулу (51) можно преобразовать аналогично, что еще более уменьшит общее количество операций. Приведем без вывода рекуррентные формулы вычисления коэффициентов Фурье при N = L^r:

,

, (52)

,

при этом .

Число вспомогательных коэффициентов k -го ранга равно N, т.е., согласно (52), для вычисления коэффициентов всех рангов требуется 2 NLr операций. Учитывая, что , можно найти оптимальное число сомножителей , при котором общее число операций 2 NLr минимально. Приравнивая производную по r выражения к нулю, находим . Таким образом, минимум числа операций 2 eN ×ln N реализуется, когда и при . Однако, в рамках процедуры быстрого преобразовании Фурье принято считать, что N = 2 ^r и L = 2. Данный случай мало отличается от оптимального, т.к. требуемое в этом случае число операций 4 N ×log₂ N мало отличается от оптимального числа операций» 3,768 N ×log₂ N.

На отрезке [0, a ] с помощью процедуры (52) быстрого преобразования Фурье численно решим задачу , правая часть которой

(53)

обеспечивает аналитическое решение вида:

. (54)

На листинге_№7 приведен код программы решения задачи (44) методом быстрого преобразования Фурье по схеме (52) и сравнение полученного численного решения с аналитическим решением (54).

Листинг_№7

%Программа решения уравнения (44) методом

%быстрого преобразования Фурье по схеме (52)

%с правой частью (53) в уравнении u''-mu u= -f

function fourier

global a mu h w

a=1; mu=1;

%Определяем возможные степени двойки числа

%узлов сетки N=2^r

r=[3 4 5 6 7 8];

%Организуем цикл расчетов с разными N=2^r

for i=1:length(r)

%Определяем параметры сетки и параметр w

N=2^r(i); h=a/N; x=0:h:a;

w=cos((2*pi)/N)+sqrt(-1)*sin((2*pi)/N);

%Вычисляем bp согласно схеме быстрого

%преобразования Фурье в (52)

for p=0:(N-1)

bp=0.5*(b(r(i),0,mod(p,2^(r(i)-1)))+...

b(r(i),1,mod(p,2^(r(i)-1)))*w^(-p));

afrr(p+1)=bp/((4/h^2)*sin((pi*p)/N)^2+mu);

end

%Находим решение согласно (45)

for n=0:N

s=0;

for q=0:(N-1)

s=s+afrr(q+1)*w^(n*q);

end

y(n+1)=s;

end

%Находим отклонение численного решения от

%аналитического, т.е. ошибку численного решения

%в норме C

for n=0:N

y(n+1)=abs(y(n+1)-x(n+1)^2*(a-x(n+1))^2);

end

error(i)=max(max(y));

Nr(i)=N;

end

%Рисуем график зависимости ошибки численного решения

%от числа шагов сетки в логарифмическом масштабе по

%осям абсцисс и ординат

loglog(Nr,error);

%Определяем функцию правой части

function y=f(x)

global a mu

y=-2*a^2+12*a*x+(mu*a^2-12)*x^2-2*mu*a*x^3+mu*x^4;

%Определяем рекурсивную функцию вычисления

%коэффициентов b(l1,l2,...,lk,pk) в быстром

%преобразовании Фурье согласно схеме (52)

function z=b(r,l,p)

global h w

k=length(l);

if k<=(r-2)

l0=[l 0]; l1=[l 1];

z=0.5*(b(r,l0,mod(p,2^(r-k)))+...

b(r,l1,mod(p,2^(r-k)))*w^(-2^k*p));

else

s=0;

for i=1:k

s=s+l(i)*2^(i-1);

end

z=0.5*(f(s*h)+f((s+2^(r-1))*h)*w^(-2^(r-1)*p));

end

На рис.7 приведен итоговый график работы кода программы листинга_№7. Из графика на рис.7 видно, что с экспоненциальный рост числа узлов сетки приводит к экспоненциальному уменьшению ошибки численного решения с помощью метода быстрого преобразования Фурье согласно схеме (52). Ошибка оценивалась по отношению к аналитическому решению (54) согласно формуле: error = .

Метод быстрого преобразования Фурье легко обобщается на многомерный случай. Для примера рассмотрим уравнение

,

для которого определена первая краевая задача в прямоугольной области. Введем равномерную сетку , n = 0,1,…, N, , m = 0,1,…, M и определим на ней разностную схему

(55)

Решение уравнения (55) будем искать в виде разложения Фурье:

.

Рис.7. Зависимость ошибки численного решения задачи (44) с
помощью быстрого преобразования Фурье (52) от N = 2 ^r

Аналогично формулам одномерного случая (46), (47) имеем

, (46¢)

где

. (47¢)

Представим (47¢) в следующем виде:

. (56)

Обе суммы в (56) имеют вид, аналогичный сумме (47). По этой причине, если N и M разлагаются на множители, то каждую из сумм можно вычислить согласно рекуррентным формулам типа (52). Если и , то число операций на каждый узел сетки есть O (r ₁ L ₁ + r ₂ L ₂) = O (log(NM)), что аналогично одномерному случаю.