9

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

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

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 ответ 1

5

К сожалению не специалист именно в теплопроводности, но могу отметить несколько точно проблемных моментов:

  1. В настоящий момент шаг по сетке времени не связан с шагом по сетке в пространстве. В явных схемах - это ведет к нестабильности. Нужно соблюдать критерий CFL. В вашем случае, если я правильно помню термодинамику, Δt < CFL * χ * (Δx)². То есть Δt < (Δx)² (χ - radiative diffusion в вашем уравнении = 1). В вашей программе это явно не соблюдается ни по конкретным значениям, ни по алгоритму. Решение: привяжите K и N друг к другу.

  2. По вашему графику очень сложно понять что происходит. Если начертить каждый шаг времени на отдельном графике - то видно, что численное решение расходится. В частности, это видно и на вашем чертеже (e121). Объясняется как минимум пунктом 1, и, возможно, еще какими-то проблемами в программе. Кроме того, решения на промежуточных шагах выглядят ну очень подозрительно.

  3. Посмотрите значения x_i и t_j (print(x_i)) - сейчас это массивы от 0 до 9.9 с шагом 0.1. Судя по условию, перепутан шаг, количество шагов и правая граница.

Теперь по проверке решений. Формально, такой код можно проверить через Method of Manufactured Solutions, но для вашего задания это, конечно, слишком.

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

  • выполняются ли начальные и конечные условия
  • анализ сходимости (убавили шаг по времени - стало лучше? убавили шаг дискретизации в пространстве - стало лучше?)
  • не появляется ли энергия из ниоткуда (для замкнутых систем)?
  • совпадает ли решение с решением из коммерческого\open-source симулятора?
  • если решить уравнение другим численным методом (например, методом конечных элементов), получим ли мы примерно тоже самое?
7
  • Не могу разобраться, что именно перепутано в пункте 3? Можно поподробнее? У меня действительно не соблюден критерий устойчивости явной схемы t/h^2<1/2. Это исправлю. Спасибо!
    – jokop
    4 ноя 2019 в 12:39
  • 1
    @jokop выведите значения массивов x_i и t_j. Они будут от 0 до 9.9. Это противоречит условию, если я правильно понимаю. Вроде как они должны быть от 0 до 1 c определенным шагом. 4 ноя 2019 в 13:03
  • Изменила K, теперь K = 25 и x_i = np.arange(0, 1, h) t_j = np.arange(0, 1, t) теперь график выглядит так: ibb.co/Qmgp555
    – jokop
    4 ноя 2019 в 13:50
  • 1
    все равно нестабильно. K=25 - маловато при N=10. Минимум K=100 (смотрите CFL), а лучше чуть больше. 4 ноя 2019 в 13:58
  • 1
    @jokop ну как минимум- он не взорвался до e121. Значит кое-какие проблемы мы уже победили :) Остальное - нужно смотреть более внимательно, с обновленным кодом и, возможно, референсным решением. 4 ноя 2019 в 14:23

Ваш ответ

By clicking “Отправить ответ”, you agree to our terms of service and acknowledge you have read our privacy policy.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.