Численное решение задачи теплопроводности
📁 python, вопрос, ошибка, пример, ии
Необходимо решить задачу теплопроводности на отрезке
При решении использовала явную схему
N = 10 # максимальное число шагов по хK = 10 # максимальное число шагов по tl = 1 # значение х на правой границеh = l / N # шаг сетки по хT = 1 # максимальное значение времени t на правой границеt = T / K # шаг сетки по времени# зададим сетку x_i = np.arange(0, N, h) # значения в узлах по хt_j = np.arange(0, K, t) # значение в узлах по tr_j = len(t_j) # количество узлов по tr_i = len(x_i) # количество узлов по xw_h_t = np.zeros([r_i, r_j]) # итоговая сетка размером x_i*t_j# зададим значение функции входящей в начальное уравнениеx = 0def f(x): retu np.sin(x)# граничные условияux_0 = 1 # граничное условие на левом конце при x=0ut_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 # максимальное число шагов по tl = 1 # значение х на правой границеh = l / N # шаг сетки по хT = 1 # максимальное значение времени t на правой границеt = T / K # шаг сетки по времени# зададим сетку x_i = np.arange(0, 1, h) # значения в узлах по хt_j = np.arange(0, 1, t) # значение в узлах по tr_j = len(t_j) # количество узлов по tr_i = len(x_i) # количество узлов по xw_h_t = np.zeros([r_j, r_i]) # итоговая сетка размером x_i*t_j# зададим значение функции входящей в начальное уравнениеx = 0def f(x): retu np.sin(x)# граничные условияux_0 = 1 # граничное условие на левом конце при x=0ut_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 К сожалению не специалист именно в теплопроводности, но могу отметить несколько точно проблемных моментов:
В настоящий момент шаг по сетке времени не связан с шагом по сетке в пространстве. В явных схемах - это ведет к нестабильности. Нужно соблюдать . В вашем случае, если я правильно помню термодинамику,
Δt < CFL * χ * (Δx)². То естьΔt < (Δx)²(χ - radiative diffusion в вашем уравнении = 1). В вашей программе это явно не соблюдается ни по конкретным значениям, ни по алгоритму. Решение: привяжитеKиNдруг к другу.По вашему графику очень сложно понять что происходит. Если начертить каждый шаг времени на отдельном графике - то видно, что численное решение расходится. В частности, это видно и на вашем чертеже (
e121). Объясняется как минимум пунктом 1, и, возможно, еще какими-то проблемами в программе. Кроме того, решения на промежуточных шагах выглядят ну очень подозрительно.Посмотрите значения
x_iиt_j(print(x_i)) - сейчас это массивы от0до9.9с шагом0.1. Судя по условию, перепутан шаг, количество шагов и правая граница.
Теперь по проверке решений. Формально, такой код можно проверить через , но для вашего задания это, конечно, слишком.
При обучении численным методам чаще всего используются задачи для которых известно аналитическое решение и полученный результат сравнивается уже с ним. Еще стоит проверять:
- выполняются ли начальные и конечные условия
- анализ сходимости (убавили шаг по времени - стало лучше? убавили шаг дискретизации в пространстве - стало лучше?)
- не появляется ли энергия из ниоткуда (для замкнутых систем)?
- совпадает ли решение с решением из коммерческого\open-source симулятора?
- если решить уравнение другим численным методом (например, методом конечных элементов), получим ли мы примерно тоже самое?
Войдите чтобы оставить ответ