Для начала рассмотрим простую одномерную задачу распространения тепла в стержне. Уравнение, описывающее распространение тепла при некотором начальном распределении температуры по стержню:
Такое уравнение решается аналитически методом разделения переменных, но нас интересует как это можно сделать численно. Прежде всего нужно определиться, как считать вторую пространственную производную по х. Проще всего это делается каким-нибудь разностным методом, например:
Но мы поступим иначе. Распределение температуры есть функция координаты и времени, и в каждый момент времени эта функция может быть представлена в виде суммы ряда Фурье, который в численном виде обрезается на n-ом члене:
Где u^«с крышечкой» — это коэффициенты разложения ряда Фурье. Подставим выражение для ряда в уравнение переноса тепла:
Получаем уравнение для коэффициентов Фурье, в котором отсутствует производная по координате! Теперь это обыкновенное дифференциальное уравнение, а не в частных производных, которое можно решить простым разностным методом. Уже легче, теперь остается найти коэффициенты разложения и в этом нам очень поможет быстрое преобразование Фурье (дальше FFT).
Логика здесь следующая:
1) в начальный момент времени дана функция координаты, описывающая распределение температуры по стержню;
2) разбиваем стержень на сетку из n точек;
3) находим комплексные коэффициенты Фурье с помощью алгоритма FFT, обозначим операцию как F(u);
4) умножаем полученные коэффиценты на -|k|2, получаем Фурье-образ второй производной. Аналогично можно получить Фурье-образ производной более высоких порядков p, достаточно умножить на (ik)p;
5) делаем обратное преобразование Фурье F-1(u), с помощью алгоритма IFFT, получаем значения второй производной в точках на сетке;
6) делаем шаг по времени, уже обычной разностной, явной или неявной, схемой;
7) повторяем.
Рассмотрим теперь как это работает в программе для Matlab/Octave. В качестве начального распределения температуры возьмем гладкую функцию u0=2+sin(x)+sin(2x), стержень длинной 2π разобьем на 50 точек, с шагом по времени h=0.1, граничные условия периодичные (кольцо).
Код для одномерного уравнения переноса тепла
clear all;n = 50; % number of points dx = 2*pi/n; % space step x = 0:dx:2*pi-dx; % grid h = 0.1; % temporal step times = 10; % number of iterations in time k = fftshift(-n/2:1:n/2-1); % wave numbers k2 = k.*k; u0 = 2 + sin(x) + sin(2*x); % initial conditions u = zeros(times,n); % stores results u(1,:) = u0; uf = fft(u0); % Fourier coefficients of initial function for i=2:times uf = uf.*(1-h*k2); % next time step in Fourier space u(i,:) = real(ifft(uf)); % IFFT to physical space end [X,T] = meshgrid(x,0:h:times*h-h);waterfall(X,T,u)
Стоит отметить особенность алгоритма FFT в Matlab, связанную с тем, что полученные коэффициенты разложения на выходе d=fft(u) идут не по порядку, а смещены, первая половина на месте второй и наоборот. Cначала идут коэффициенты с номерами от 0 до n/2-1, потом с номерами от -n/2 до -1. С этим были проблемы…
Полученное решение можно видеть на графике в виде «водопада» линий распределения температуры по х для каждого момента времени t. Видно, что решение испытывает численную неустойчивость, связано это с невыполнением критерия Куранта. Избавиться от неустойчивости можно уменьшив шаг по времени, либо применяя более продвинутую неявную схему, например Кранка-Николсона [3]
Заключение
В данной работе рассматриваются методы решения одной из базовых задач, рассматриваемых в курсе математической физики, использующих уравнения распространения тепла в стержне. В применении к уравнениям в частных производных схема разделения переменных приводит к нахождению решения в виде ряда или интеграла Фурье. В этом случае метод также называют методом Фурье (в честь Жана Батиста Фурье, построившего решения уравнения теплопроводности в виде тригонометрических рядов) и методом стоячих волн.
В работе рассматриваются как аналитический метод решения, где результат выводится при помощи математических преобразований, так и численный метод, при котором результат соответствует действительному с заданной точностью, но который требует много рутинных вычислений и поэтому выполним только при помощи вычислительной техники (ЭВМ). В качестве числового расчета на ЭВМ в работе рассматривается так называемый спектральный метод, использующий программу Mathlab.
Основная суть спектрального метода, это замена исходных дифференциальных уравнений в частных производных на обыкновенные диффуры для коэффициентов разложения искомых функций по некоторому базису. Базисом могут быть синусы-косинусы, комплексные экспоненты, ортогональные полиномы, если требует геометрия — цилиндрические или сферические функции. Найденные коэффициенты в каждый момент времени позволяют восстановить искомое решение, а алгоритм FFT позволяет делать это быстро. Спектральный метод эффективен в отношении памяти, поэтому широко используется в геофизике, моделировании климата и метеорологии.
Список источников
- Лободенко Е.И., Курс лекции «Уравнения математической физики»
- https://x-hunt.ru/zadasha/algebra37.html «Введение в математический анализ». Информационный портал.
- https://habrahabr.ru/post/267401/ Информационный портал для IT специалистов. Статья «Спектральный метод на основе простых задач матфизики»
- Сабитов, К.Б. Уравнения математической физики // Москва: Физматлит. – 2013