Численное решение задачи теплопроводности

#python #численные_методы

Необходимо решить задачу теплопроводности на отрезке


При решении использовала явную схему

N = 10 # максимальное число шагов по х
K = 10 # максимальное число шагов по t
l = 1 # значение х на правой границе
h = l / N # шаг сетки по х
T = 1 # максимальное значение времени t на правой границе
t = T / K # шаг сетки по времени

# зададим сетку 
x_i = np.arange(0, N, h) # значения в узлах по х
t_j = np.arange(0, K, t) # значение в узлах по t
r_j = len(t_j) # количество узлов по t
r_i = len(x_i) # количество узлов по x
w_h_t = np.zeros([r_i, r_j]) # итоговая сетка размером x_i*t_j

# зададим значение функции входящей в начальное уравнение
x = 0
def f(x):
    return np.sin(x)

# граничные условия
ux_0 = 1 # граничное условие на левом конце при x=0
ut_0 = np.cos(x_i) # граничное условие при t=0

# найдем значения на нулевом слое при t=0 ut_0 = np.cos(x_i)
w_h_t[0] = np.cos(x_i)

# найдем значения w_h_t на первом и последующих слоях
const = t / (h**2) 
for j in range(1, len(x_i)-1):
    for i in range(len(w_h_t[j])-1):        
        w_h_t[j+1, i] = w_h_t[j, i] + const * (w_h_t[j,i] - 2*w_h_t[j,i] + w_h_t[j,
i-1]) + t*f(x_i[j])
        w_h_t[j+1, 0] = 1
        w_h_t[j+1, len(w_h_t[i])-1] = w_h_t[j+1, len(w_h_t[i])-2] + h * t_j[j+1]


plot_ = np.arange(0,len(w_h_t)-1,1)
for y in plot_:
    plt.plot(x_i, w_h_t[y])


В результате получила вот такой график зависимости Х от рассчитанного значения функции
в узлах сетки 



Мне не понятно насколько неверно мое решение. Возникли проблемы с поиском частного
решения, wolfram выдал u(x) = 1.54x+1+sinx, что мне кажется не верным, а самой решить
не получилось. В учебнике Филиппова похожих примеров не нашлось, в сети ничего достаточно
подробного, чтобы разобраться в решении не нашла. Подскажите где можно найти как решать
аналитически такое уравнение и насколько неверно мое решение? В чем ошибка? И как вообще
проверяют на корректность решения таких задач, кроме как сравнения с аналитическим
решением?



В общем и целом разобралась. Решила дополнить свой вопрос отредактированным решением,
возможно, кому-то пригодится.
Численное решение правильно найти так и не удалось, но находится оно методом Фурье(разделение
переменных). 
Итоговый график:
 

Рабочий код:

N = 10 # максимальное число шагов по х
K = 500 # максимальное число шагов по t
l = 1 # значение х на правой границе
h = l / N # шаг сетки по х
T = 1 # максимальное значение времени t на правой границе
t = T / K # шаг сетки по времени


# зададим сетку 
x_i = np.arange(0, 1, h) # значения в узлах по х
t_j = np.arange(0, 1, t) # значение в узлах по t
r_j = len(t_j) # количество узлов по t
r_i = len(x_i) # количество узлов по x
w_h_t = np.zeros([r_j, r_i]) # итоговая сетка размером x_i*t_j


# зададим значение функции входящей в начальное уравнение
x = 0
def f(x):
    return np.sin(x)

# граничные условия
ux_0 = 1 # граничное условие на левом конце при x=0
ut_0 = np.cos(x_i) # граничное условие при t=0

# найдем значения на нулевом слое при t=0 ut_0 = np.cos(x_i)
w_h_t[0] = np.cos(x_i)

# найдем значения w_h_t на первом и последующих слоях
const = t / (h**2) 
for j in range(len(w_h_t) - 1):
    for i in range(len(w_h_t[j]) - 1):        
        w_h_t[j + 1, i] = w_h_t[j, i] + const* (w_h_t[j, i+1] - 2 * w_h_t[j, i] +
w_h_t[j, i - 1]) + t*f(x_i[i])
        w_h_t[j + 1, 0] = 1
        w_h_t[j + 1, len(w_h_t[i])-1] = w_h_t[j + 1, len(w_h_t[i])-1] + h

    

Ответы

Ответ 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, но для вашего задания это, конечно, слишком. 

При обучении численным методам чаще всего используются задачи для которых известно
аналитическое решение и полученный результат сравнивается уже с ним. Еще стоит проверять:


выполняются ли начальные и конечные условия
анализ сходимости (убавили шаг по времени - стало лучше? убавили шаг дискретизации
в пространстве - стало лучше?)
не появляется ли энергия из ниоткуда (для замкнутых систем)?
совпадает ли решение с решением из коммерческого\open-source симулятора?
если решить уравнение другим численным методом (например, методом конечных элементов),
получим ли мы примерно тоже самое?