Разработать компьютерную программу для численного решения уравнения теплопроводности
Привет, Хабр! Некоторое время назад увлекся глубоким обучением и стал потихоньку изучать tensorflow. Пока копался в tensorflow вспомнил про свою курсовую по параллельному программированию, которую делал в том году на 4 курсе университета. Задание там формулировалось так:
Линейная начально-краевая задача для двумерного уравнения теплопроводности:
Хотя правильнее было бы назвать это уравнением диффузии.
Задачу тогда требовалось решить методом конечных разностей по неявной схеме, используя MPI для распараллеливания и метод сопряженных градиентов.
Я не специалист в численных методах, пока не специалист в tensorflow, но опыт у меня уже появился. И я загорелся желанием попробовать вычислять урматы на фреймворке для глубокого обучения. Метод сопряженных градиентов реализовывать второй раз уже не интересно, зато интересно посмотреть как с вычислением справится tensorflow и какие сложности при этом возникнут. Этот пост про то, что из этого вышло.
Численный алгоритм
Разностная схема:
Чтобы проще было расписывать, введем операторы:
Явная разностная схема:
В случае явной разностной схемы для вычисления используются значения функции в предыдущий момент времени и не требуется решать уравнение на значения . Однако такая схема менее точная и требует значительно меньший шаг по времени.
Неявная разностная схема:
Перенесем в левую сторону все связанное с , а в правую и домножим на :
По сути мы получили операторное уравнение над сеткой:
что, если записать значения в узлах сетки как обычный вектор, является обычной системой линейных уравнений (). Значения в предыдущий момент времени константы, так как уже рассчитаны.
Для удобства представим оператор как разность двух операторов:
Заменив на нашу оценку , запишем функционал ошибки:
где — ошибка в узлах сетки.
Будем итерационно минимизировать функционал ошибки, используя градиент.
В итоге задача свелась к перемножению тензоров и градиентному спуску, а это именно то, для чего tensorflow и был задуман.
Реализация на tensorflow
Кратко о tensorflow
В tensorflow сначала строится граф вычислений. Ресурсы под граф выделяются внутри tf.Session. Узлы графа — это операции над данными. Ячейками для входных данных в граф служат tf.placeholder. Чтобы выполнить граф, надо у объекта сессии запустить метод run, передав в него интересующую операцию и входные данные для плейсхолдеров. Метод run вернет результат выполнения операции, а также может изменить значения внутри tf.Variable в рамках сессии.
tensorflow сам умеет строить графы операций, реализующие backpropagation градиента, при условии, что в оригинальном графе присутствуют только операции, для которых реализован градиент (пока не у всех).
Сначала код инициализации. Здесь производим все предварительные операции и считаем все, что можно посчитать заранее.
и следует брать такими, чтобы было небольшим, желательно, хотя бы < 1, особенно при использовании «негладких» функций.
Метод который строит граф уравнения:
По-хорошему надо было считать значения функции на краях заданными и оптимизировать значения функции только во внутренней области, но с этим возникли проблемы. Способа сделать оптимизируемым только часть тензора не нашлось, и у операции присвоения значения срезу тензора не написан градиент (на момент написания поста). Можно было бы попробовать хитро повозиться на краях или написать свой оптимизатор. Но и просто добавление разности на краях значений функции и краевых условий в функционал ошибки хорошо работает.
Стоит отметить, что метод с адаптивным моментом показал себя наилучшим образом, пусть функционал ошибки и квадратичный.
Вычисление функции: в каждый момент времени делаем несколько оптимизационных итераций, пока не превысим maxiter или ошибка не станет меньше eps, сохраняем и переходим к следующему моменту.
Вернемся к нашей чашечке кофе и разработаем компьютерную программу для численного решения уравнения теплопроводности Ньютона. Чтобы определить, разумна ли наша модель понижения температуры, приведем экспериментальные данные для остывания настоящей чашки кофе (табл. 2.2.)
Остывание чашки кофе, помещенной
в керамический стакан.
Время, мин | Т, град | Время, мин | Т, град |
0 | 83.0 | 8.0 | 64.7 |
1.0 | 77.7 | 9.0 | 63.4 |
2.0 | 75.1 | 10.0 | 62.1 |
3.0 | 73.0 | 11.0 | 61.0 |
4.0 | 71.1 | 12.0 | 59.9 |
5.0 | 69.4 | 13.0 | 58.7 |
6.0 | 67.8 | 13.0 | 57.8 |
7.0 | 66.4 | 15.0 | 56.6 |
Температура регистрировалась с точностью 0.1 град.
Температура окружающего воздуха 22 град.
Структура программы Сool, в которой реализован алгоритм Эйлера численного решения задачи ньютоновской теплопроводности, аналогична структуре программы Example. Единственное различие состоит в том, что вывод на экран осуществляется подпрограммой Output, которая вызывается из подпрограммы Euler. Модульный подход к программированию будет полезен при выводе графиков.
program Cool;
t, tmpr, room_tmpr, r, dt : real;
procedure Iinitial ( var t, tmpr, room_tmpr, r, dt : real;
var ncalc : integer);
procedure Output(t,tmpr : real);
procedure Euler (var t, tmpr, room_tmpr, r, dt : real;
var ncalc : integer);
for icalc:=1 to ncalc do
change := - r *(tmpr - room_tmpr);
tmpr := tmpr + change* dt;
Output (t,tmpr)
Initial(t, tmpr, room_tmpr, r, dt, ncalc);
Euler(t, tmpr, room_tmpr, r, dt, ncalc) ;
Задача 2.1. Программа для решения задачи об остывании кофе
а). После того набора текста программы и устранения всех синтаксических ошибок, необходимо убедиться в правильной реализации в программе требуемого алгоритма. Самое простое, что можно сделать - это сравнить численные результаты с предельными случаями, для которых имеется аналитическое решение, или вычисления можно проделать вручную. Используйте метод Эйлера и калькулятор для численного решения задачи теплопроводности Ньютона с теми же параметрами, что и в программе Сool. Сравните свои расчеты с тем, что получилось по программе Cool, и проверьте свою программу в этом случае.
б). Модифицируйте программу Coolтак, чтобы значения параметров r, dt, tmax можно было вводить с клавиатуры, а также, чтобы перед таблицей, содержащей время и температуру, печатался заголовок.
Поскольку приведенная задача программирования может оказаться для вас новой, ниже приводится “решение”. Заметим, что в модифицированной подпрограмме Euler производится ncalc итерации между обращениями к подпрограмме Оutput. Конечно же, наше решение не единственное и должно рассматриваться только в качестве одного из вариантов.
Как указывалось ранее, точное решение дифференциального уравнения теплопроводности (1.22) удается получить лишь для немногих частных случаев. Поэтому на практике применяются различные приближенные модели решений. Численные методы описания тепловых процессов находят все более широкое применение в связи с развитием вычислительной техники.
Наибольшее распространение в решениях получили модели, основанные на методе сеток. Сущность метода заключается в том, что исследуемая непрерывная функция заменяется совокупностью приближенных значений, рассчитанных в некоторых точках-узлах. Совокупность узлов, соединенных определенным образом, образует сетку. Сетка является дискретной моделью исследуемой непрерывной функции.
Применение метода сеток позволяет свести решение задачи к решению системы алгебраических уравнений, что легко осуществляется на ЭВМ. При решении теплофизических задач наибольшее применение получили следующие виды метода сеток – метод конечных разностей (МКР), метод конечных элементов (МКЭ), метод граничных элементов (МГЭ).
Метод конечных разностей (МКР). Предусматривает разбиение области определения функции на отдельные равные элементы. Рассмотрим, например, распределение температуры вдоль стержня малого сечения (рис.3.4.)
Рис. 3.4 – Использование метода конечных разностей для расчета температурного поля в стержне
Температурное поле в таком стержне будет одномерным (см.процесс иглофрезерования).
Разделим стержень на n равных элементов размером hх.
Время процесса также разделим на равные промежутки Δτ с нумерацией 1, 2, 3. (k + 1), k, (k – 1). m, т.е. m промежутков.
Можно доказать, что
Формула (3.9) выражает дифференциальное уравнение теплопроводности, представленное в конечно-разностной форме. Как видно, это уравнение алгебраическое. Чтобы определить температуру в любой точке тела в данный момент времени (i, k + 1) достаточно знать температуру соседних точек (i+1, i – 1) в предыдущий момент времени (k). Используя граничные начальные условия (они заданы), получаем систему связанных друг с другом уравнений. Таких уравнений будет nm с таким же числом неизвестных. Решая эту систему можно рассчитать температуру в каждой точке твердого тела в любой момент времени. Результат получается тем точнее, чем меньше будут размер элементов hx и промежуток времени Δτ.
Метод конечных элементов (МКЭ). В сложных случаях теплообмена разбиение твердого тела на одинаковые по размеру элементы, а времени – на одинаковые промежутки вызывает большой объем вычислительной работы. В любой конкретной задаче разные участки твердых тел представляют для практики различный интерес. Некоторые объемы твердых тел изучаются с большей степенью детализации (например, температура контактных поверхностей режущих инструментов), чем другие. Метод конечных элементов позволяет осуществлять различную детализацию решения в различных областях изучаемого объекта, причем могут быть использованы элементарные объемы различные не только по величине, но и по конфигурации. Однако решение задачи в этом случае выполняется с использованием вариационного исчисления. При использовании МКЭ число уравнений и число неизвестных может быть значительно меньше, чем при методе конечных разностей, однако математический аппарат существенно сложнее.
Метод граничных элементов (МГЭ)
Существенное отличие МГЭ от МКР и МКЭ состоит в том, что на конечные элементы разбивают только граничные поверхности тела, а не весь его объем. Реализация этого метода приводит к значительным математическим трудностям, поэтому этот метод при анализе тепловых процессов пока применяют значительно редко.
В настоящее время разработано множество методов решения дифференциального уравнения теплопроводности: аналитические, численные и методы математического моделирования.
Решение дифференциального уравнения теплопроводности при использовании аналитических методов получается в виде формулы, а для численных методов − в виде массива чисел.
Достоинство аналитических методов: обозримость получаемых результатов, возможность полного анализа полученного решения. Недостаток аналитических методов: возможность получения решения, как правило, для областей канонической формы. Для декартовой системы координат область канонической формы ограничивается линиями, параллельными осям OX и OY.
Достоинство численных методов: возможность решать сложные задачи, недоступные для аналитических методов. Массив чисел, полученных численными методами, обычно представляется в табличной форме или аппроксимируется формулой для последующего анализа.
Аналитические и численные методы следует рассматривать как дополняющие друг друга, а не как конкурирующие.
Для решения задач теплофизики резания наибольшее распространение получили следующие аналитические методы: метод точечных источников и метод быстродвижущихся источников тепла. Рассмотрим их по порядку.
Решением дифференциального уравнения теплопроводности (1.15) является функция , которая при подстановке в уравнение (1.15) обращает его в тождество и, кроме того, удовлетворяет краевым условиям.
Пусть в бесконечном стержне с нулевой начальной температурой, (рис. 1.2) в начальный момент времени в точке x=x 1 вспыхнул и мгновенно погас источник тепла, выделивший количество тепла Q.
Граничные условия могут быть записаны в виде
Это указывает на отсутствие теплообмена между стержнем и окружающей средой на бесконечности.
Решение уравнения теплопроводности для мгновенного точечного источника тепла предложено У.Кельвином и имеет вид [1]
Непосредственной проверкой можно убедиться, что функция (1.20) удовлетворяет уравнению теплопроводности (1.15) и граничным условиям (1.19). Из (1.20) следует, что функция имеет максимум в точке x=x 1 и что количество тепла Q в любой момент времени остаётся неизменным и равным Cv·В, а также, что величина B представляет собой площадь, ограниченную функцией и осью x.
Функцию называют фундаментальным решением уравнения теплопроводности. С помощью этого уравнения можно сконструировать решение
уравнения теплопроводности для других краевых условий: для этого любой процесс распространения тепла в твёрдом теле теплопроводностью представляется как совокупность процессов выравнивания температуры от множества элементарных источников тепла, распределённых как в пространстве, так и во времени. Результат получается суммированием (интегрированием) элементарных решений.
Рассмотрим метод быстродвижущих источников тепла. Быстродвижущими называются такие источники, скорость перемещения которых больше скорости распространения тепла. К быстродвижущимся источникам относят такие, для которых критерий Пекле: , рис.1.3.
Здесь полуплоскость , имеющая температуру , движется с постоянной скоростью V относительно оси y. На границе полуплоскости x=0 расположен равномерно распределённый теплоисточник, q=const, длина которого равна C. Суть метода быстродвижущихся источников тепла заключается в следующем: в среде, движущихся с большой скоростью V вдоль оси y относительно теплоисточника можно пренебречь переносом тепла за счёт теплопроводности вдоль этой оси.
Из закона Фурье следует, что а значит . Таким образом, задача значительно упрощается, т.к. дифференциальное уравнение теплопроводности для двумерного поля (1.16) сводится к дифференциальному уравнению теплопроводности для одномерного поля (1.17) в полубесконечном стержне (заштрихованном) (рис.1.3), движущемся вдоль оси y относительно источника тепла q=const.
Необходимо решить задачу теплопроводности на отрезке
При решении использовала явную схему
В результате получила вот такой график зависимости Х от рассчитанного значения функции в узлах сетки
Мне не понятно насколько неверно мое решение. Возникли проблемы с поиском частного решения, wolfram выдал u(x) = 1.54x+1+sinx , что мне кажется не верным, а самой решить не получилось. В учебнике Филиппова похожих примеров не нашлось, в сети ничего достаточно подробного, чтобы разобраться в решении не нашла. Подскажите где можно найти как решать аналитически такое уравнение и насколько неверно мое решение? В чем ошибка? И как вообще проверяют на корректность решения таких задач, кроме как сравнения с аналитическим решением?
В общем и целом разобралась. Решила дополнить свой вопрос отредактированным решением, возможно, кому-то пригодится. Численное решение правильно найти так и не удалось, но находится оно методом Фурье(разделение переменных). Итоговый график:
1 ответ 1
К сожалению не специалист именно в теплопроводности, но могу отметить несколько точно проблемных моментов:
В настоящий момент шаг по сетке времени не связан с шагом по сетке в пространстве. В явных схемах - это ведет к нестабильности. Нужно соблюдать критерий CFL. В вашем случае, если я правильно помню термодинамику, Δt < CFL * χ * (Δx)² . То есть Δt < (Δx)² (χ - radiative diffusion в вашем уравнении = 1). В вашей программе это явно не соблюдается ни по конкретным значениям, ни по алгоритму. Решение: привяжите K и N друг к другу.
По вашему графику очень сложно понять что происходит. Если начертить каждый шаг времени на отдельном графике - то видно, что численное решение расходится. В частности, это видно и на вашем чертеже ( e121 ). Объясняется как минимум пунктом 1, и, возможно, еще какими-то проблемами в программе. Кроме того, решения на промежуточных шагах выглядят ну очень подозрительно.
Посмотрите значения x_i и t_j ( print(x_i) ) - сейчас это массивы от 0 до 9.9 с шагом 0.1 . Судя по условию, перепутан шаг, количество шагов и правая граница.
Теперь по проверке решений. Формально, такой код можно проверить через Method of Manufactured Solutions, но для вашего задания это, конечно, слишком.
При обучении численным методам чаще всего используются задачи для которых известно аналитическое решение и полученный результат сравнивается уже с ним. Еще стоит проверять:
Читайте также: