Я решаю уравнение теплопроводности для металлического стержня, поскольку один конец поддерживается при 100 ° C, а другой при 0 ° C как
Если изменение граничного условия Неймана как одного конца изолировано (а не флюсом),
то как расчетный термин
Типичный подход к граничному условию Неймана состоит в том, чтобы представить «точку-призрак» на один шаг за пределами области и вычислить значение для нее с использованием граничного условия; затем выполните нормально (используя PDE) для точек, находящихся внутри сетки, включая границу Неймана.
Точка призрака позволяет использовать симметричную конечную разностную аппроксимацию производной на границе, т.е. (T[n, j+1] — T[n, j-1])/(2*dy) если y — пространственная переменная. Несимметричная аппроксимация (T[n, j] — T[n, j-1])/dy , которая не содержит точки призрака, намного менее точна: введенная ошибка на порядок хуже ошибки участвующих в дискретизации самого PDE.
Итак, когда j является максимальным возможным индексом для T, граничное условие говорит, что » T[n, j+1] » следует понимать как T[n, j-1] и это то, что делается ниже.
Видео:Python - численное решение дифференциального уравнения 1го порядка и вывод графикаСкачать
Уравнение теплопроводности в tensorflow
Привет, Хабр! Некоторое время назад увлекся глубоким обучением и стал потихоньку изучать tensorflow. Пока копался в tensorflow вспомнил про свою курсовую по параллельному программированию, которую делал в том году на 4 курсе университета. Задание там формулировалось так:
Линейная начально-краевая задача для двумерного уравнения теплопроводности:
Хотя правильнее было бы назвать это уравнением диффузии.
Задачу тогда требовалось решить методом конечных разностей по неявной схеме, используя MPI для распараллеливания и метод сопряженных градиентов.
Я не специалист в численных методах, пока не специалист в tensorflow, но опыт у меня уже появился. И я загорелся желанием попробовать вычислять урматы на фреймворке для глубокого обучения. Метод сопряженных градиентов реализовывать второй раз уже не интересно, зато интересно посмотреть как с вычислением справится tensorflow и какие сложности при этом возникнут. Этот пост про то, что из этого вышло.
Видео:#5. Математические функции и работа с модулем math | Python для начинающихСкачать
Численный алгоритм
Разностная схема:
Чтобы проще было расписывать, введем операторы:
Явная разностная схема:
В случае явной разностной схемы для вычисления используются значения функции в предыдущий момент времени и не требуется решать уравнение на значения . Однако такая схема менее точная и требует значительно меньший шаг по времени.
Неявная разностная схема:
Перенесем в левую сторону все связанное с , а в правую и домножим на :
По сути мы получили операторное уравнение над сеткой:
что, если записать значения в узлах сетки как обычный вектор, является обычной системой линейных уравнений (). Значения в предыдущий момент времени константы, так как уже рассчитаны.
Для удобства представим оператор как разность двух операторов:
Заменив на нашу оценку , запишем функционал ошибки:
где — ошибка в узлах сетки.
Будем итерационно минимизировать функционал ошибки, используя градиент.
В итоге задача свелась к перемножению тензоров и градиентному спуску, а это именно то, для чего tensorflow и был задуман.
Видео:Решение уравнения теплопроводности методом конечных разностейСкачать
Реализация на tensorflow
Кратко о tensorflow
В tensorflow сначала строится граф вычислений. Ресурсы под граф выделяются внутри tf.Session. Узлы графа — это операции над данными. Ячейками для входных данных в граф служат tf.placeholder. Чтобы выполнить граф, надо у объекта сессии запустить метод run, передав в него интересующую операцию и входные данные для плейсхолдеров. Метод run вернет результат выполнения операции, а также может изменить значения внутри tf.Variable в рамках сессии.
tensorflow сам умеет строить графы операций, реализующие backpropagation градиента, при условии, что в оригинальном графе присутствуют только операции, для которых реализован градиент (пока не у всех).
Сначала код инициализации. Здесь производим все предварительные операции и считаем все, что можно посчитать заранее.
По-хорошему надо было считать значения функции на краях заданными и оптимизировать значения функции только во внутренней области, но с этим возникли проблемы. Способа сделать оптимизируемым только часть тензора не нашлось, и у операции присвоения значения срезу тензора не написан градиент (на момент написания поста). Можно было бы попробовать хитро повозиться на краях или написать свой оптимизатор. Но и просто добавление разности на краях значений функции и краевых условий в функционал ошибки хорошо работает.
Стоит отметить, что метод с адаптивным моментом показал себя наилучшим образом, пусть функционал ошибки и квадратичный.
Вычисление функции: в каждый момент времени делаем несколько оптимизационных итераций, пока не превысим maxiter или ошибка не станет меньше eps, сохраняем и переходим к следующему моменту.
Запуск:
Видео:Python для самых маленьких. Линейные уравнения. Решение задачСкачать
Результаты
Условие как и оригинальное, но без в уравнении:
Что легко правится в коде:
Разницы почти нет, потому что производные имеют большие порядки, чем сама функция.
Условие с одним нагревающимся краем:
Условие с остыванием изначально нагретой области:
Условие с включением нагрева в области:
Видео:01.02. Модель SIR. Численное решение системы дифференциальных уравнений с помощью SciPyСкачать
Рисование гифок
Функция рисования 3D-гифки:
В основной класс добавляем метод, возвращающий U в виде pandas.DataFrame
Функция рисования 2D-гифки:
Стоит отметить, что оригинальное условие без использования GPU считалось 4м 26с, а с использованием GPU 2м 11с. При больших значениях точек разрыв растет. Однако не все операции в полученном графе GPU-совместимы.
- Intel Core i7 6700HQ 2600 МГц,
- NVIDIA GeForce GTX 960M.
Посмотреть, какие операции на чем выполняются, можно с помощью следующего кода:
Это был интересный опыт. Tensorflow неплохо показал себя для этой задачи. Может быть даже такой подход получит какое-то применение — всяко приятнее писать код на питоне, чем на C/C++, а с развитием tensorflow станет еще проще.
Видео:Лекция №1.1 Явная и неявная схемы для уравнения теплопроводностиСкачать
Научный форум dxdy
Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки
| Вход Регистрация | Donate FAQ Правила Поиск |
Правила форума
Видео:Использование библиотеки SymPy для работы с системами уравнений в PythonСкачать
Уравнение теплопроводности. Продольно-поперечная схема.
 |
Последний раз редактировалось Pphantom 26.11.2017, 22:35, всего редактировалось 1 раз.
Стоит задача решить задачу дирехле для уравнения теплопроводности (в прямоугольной области). Насколько я знаю, продольно-поперечная схема является безусловно устойчивой, значит должна сходиться при любом шаге по времени или координате. Я попытался реализовать алгоритм, в итоге никакой безусловной устойчивости не получается. Ведёт себя, как явная схема: если уменьшить шаг по координатам, то, без уменьшения шага по времени, решение расходится. С чем это может быть связано? Я думал, может просто ошибка машинном округлении накапливается.
Если что, вот питоний код моей реализации. Функция принимает значения на одном временном слое и находит значения на следующем.
#
# T — шаг по времени
# hx = hy шаги по координатам
#
def solveLayer ( vv ) :
for k in range ( 1 , Nxy + 1 ) : # Первый полушаг по времени
#Заполнение коэффициентов СЛУ для неявной схемы. Крайние элементы VV всегда = 0 (гран условия)
a = [ 0 if i == 1 else — 0.5 * T / ( hx ** 2 ) for i in range ( 1 , Nxy + 1 ) ]
b = [ T / ( hx ** 2 ) + 1 ] * Nxy
c = [ 0 if i == Nxy else — 0.5 * T / ( hx ** 2 ) for i in range ( 1 , Nxy + 1 ) ]
d = [ vv [ i ] [ k ] * ( 1 — T / ( hy ** 2 ) ) + ( vv [ i ] [ k + 1 ] + vv [ i ] [ k — 1 ] ) * 0.5 * T / ( hy ** 2 ) for i in
range ( 1 , Nxy + 1 ) ]
for j in range ( 1 , Nxy ) :
t = a [ j ] / b [ j — 1 ]
a [ j ] = 0
b [ j ] — = t * c [ j — 1 ]
d [ j ] — = t * d [ j — 1 ]
vv [ Nxy ] [ k ] = d [ Nxy — 1 ] / b [ Nxy — 1 ]
for i in range ( Nxy — 1 , 0 , — 1 ) :
vv [ i ] [ k ] = 1 / b [ i — 1 ] * ( d [ i — 1 ] — c [ i — 1 ] * vv [ i + 1 ] [ k ] )
for k in range ( 1 , Nxy + 1 ) : # Второй полушаг
a = [ 0 if i == 1 else — 0.5 * T / ( hx ** 2 ) for i in range ( 1 , Nxy + 1 ) ]
b = [ T / ( hx ** 2 ) + 1 ] * Nxy
c = [ 0 if i == Nxy else — 0.5 * T / ( hx ** 2 ) for i in range ( 1 , Nxy + 1 ) ]
d = [ vv [ k ] [ i ] * ( 1 — T / ( hy ** 2 ) ) + ( vv [ k + 1 ] [ i ] + vv [ k — 1 ] [ i ] ) * 0.5 * T / ( hy ** 2 ) for i in
range ( 1 , Nxy + 1 ) ]
for j in range ( 1 , Nxy ) :
t = a [ j ] / b [ j — 1 ]
a [ j ] = 0
b [ j ] — = t * c [ j — 1 ]
d [ j ] — = t * d [ j — 1 ]
vv [ k ] [ Nxy ] = d [ Nxy — 1 ] / b [ Nxy — 1 ]
for i in range ( Nxy — 1 , 0 , — 1 ) :
vv [ k ] [ i ] = 1 / b [ i — 1 ] * ( d [ i — 1 ] — c [ i — 1 ] * vv [ k ] [ i + 1 ] )
return vv
| i | Pphantom: |
| Очень полезно пользоваться подсветкой кода, благо она есть. Поправил. |
Нет, этого не может быть. Дело в том, что схема переменных направлений — она, в общем, уточняющая по своему смыслу. Во всяком случае, через каждую пару шагов, а этого и достаточно. Поэтому погрешности округления будут, наоборот, нивелироваться. Если бы речь шла о задаче Неймана (скажем), то эффект мог бы объясняться неудачной аппроксимакцией граничных условий. Но для задачи Дир и хле в прямоугольнике никакие артефакты невозможны. Какая-то ошибка в коде. | ||
| ||
01/08/06 |
| ||||||||
где
и все эти коэффициенты положительны.
