В этой статье рассмотрим, как можно посчитать в Mathcad Express дифференциальное уравнение диффузии тепла при помощи самой простой явной разностной схемы и остановимся на свойстве устойчивости разностных схем. Речь пойдет о том, как получить вот такое решение уравнения, моделирующего остывание одномерного объекта:
На графике показано начальное распределение температуры вдоль оси Х (красная линия) и два расчетных профиля – после первого шага и после нескольких шагов по времени.
Уравнение теплопроводности
В качестве примера дифференциального уравнения в частных производных мы рассматриваем уравнение теплопроводности (или, по-другому, диффузии тепла).
Это уравнение параболического типа, содержащее первую производную по времени t и вторую по пространственной координате x. Оно описывает динамику температуры u(x,t) в присутствии источников тепла, например, при нагреве металлического стержня. Таким образом, неизвестной функцией, подлежащей определению, является функция u(x,t). Это уравнение линейно, если источник и коэффициент диффузии D либо не зависят от температуры, либо источник зависит от нее линейно, т.к. в этом случае неизвестная функция u(x,t) (и все ее производные) входит в уравнение в первой степени.
Линейное уравнение теплопроводности имеет аналитическое решение, в то время, как подавляющее большинство нелинейных уравнений в частных производных можно решить только численно. Для того, чтобы уравнения в частных производных имело единственное решение, необходимо поставить нужное количество начальных и граничных условий, т.е., в данном случае, соотношения типа u(x,0)=Init(x) (начальное условие) вместе с u(0,t)=f1(t) и u(L,t)=f2(t) (граничные условия, если уравнение решается на интервале от 0 до L).
В приведенной интерпретации одномерное уравнение диффузии тепла описывает динамику температуры некоторого удлиненного тела, например, металлического стержня:
Мы будем рассматривать случай нулевого источника тепла, а в качестве начального условия выберем некоторый реалистичный профиль температуры (стержень, сильно нагретый в центре). В качестве граничных условий зададим поддержание краев стержня при постоянной температуре. Разумеется, динамика решения описывает ожидаемую картину выравнивания температуры стержня со временем (за счет перетекания тепла от центра к периферии стержня благодаря механизму теплопроводности). Скорость передачи тепла определяется значением коэффициента диффузии D. (Заинтересовавшемуся читателю я предлагаю поэкспериментировать с прилагаемыми расчетами в Mathcad Express, выбирая различные значения коэффициента диффузии).
Разностная схема
Покажем, как реализовывать в Mathcad основной вычислительный подход к решению уравнений в частных производных – метод сеток. Будем делать это «вручную», чтобы не выйти за рамки бесплатной версии Mathcad Express.
Суть метода сеток заключается в покрытии расчетной области (x,t) сеткой из МхN точек, что определит шаги по времени и пространству ? и ? соответственно. Тем самым, определяются узлы, в которых будет осуществляться поиск решения. Затем надо заменить дифференциальное уравнения аппроксимирующим его уравнением в конечных разностях, для каждого (i,n)-го узла сетки. В случае уравнения теплопроводности достаточно просто заменить первую производную по времени и вторую по пространству их разностными аналогами (такой метод дискретизации называют методом Эйлера, а получившуюся систему разностных уравнений – разностной схемой).
Выписанная схема является явной, т.к. в ней неизвестное значение искомой функции (на n+1 шаге по времени) можно выразить через уже вычисленные ранее значения функции на предыдущем n-м шаге («слое») по времени. Более того, т.к. и на каждом слое уравнения связаны между собой рекуррентно, т.е. значение в каждом неизвестном узле на n-м слое можно найти, зная только вычисленные ранее значения с предыдущего слоя. Это позволяет организовать процесс т.н. «бегущего счета», решая соответствующие уравнения на каждом слое (исходя из известного решения на предыдущем слое).
Следующий расчет реализует описанный алгоритм разностной явной схемы бегущего счета в Mathcad Express:
Собственно, вычисленные значения матрицы U и приведены на графике, который был показан в самом начале статьи (до ката).
Устойчивость
В заключение, остановимся на свойстве устойчивости разностных схем. Параметр задачи, который характеризует отношение шагов разностной схемы по пространству и времени, называется коэффициентом (или «числом») Куранта:
Существенно, что предыдущие расчеты по явной схеме проводились при соотношении Куранта менее 1. А вот если попробовать увеличить шаг по времени в явной схеме, так, чтобы соотношение Куранта стало больше 1, то мы столкнемся с катастрофой: решение «разваливается», и вместо ожидаемого (и просчитанного при меньших соотношениях Куранта случаях) решения получается быстроосцилирующее, совершенно неправдоподобное решение.
Этот пример демонстрирует важное свойство разностных схем, называемое устойчивостью. Оказывается, многие схемы (безусловно, или при определенном соотношении параметров, как наша схема) им не обладают, и применять такие схемы для решения задач нельзя.
Более подробную информацию о разностных схемах и примерах их расчета в Mathcad читатель может найти в моей книге по вычислительной физике (разностным схемам посвящены главы 5 и 6, книга в свободном доступе). В частности, предлагаю попробовать решить разностным методом «обратное уравнение теплопроводности», которое получается заменой переменной t на -t и которое (по идее) должно описывать реконструкцию исходного профиля температуры стержня по данным наблюдения через некоторое время.
И повторюсь, что файл с расчетами, описанными в этой статье, вы найдете здесь, а дистрибутив Mathcad Express можно скачать по ссылке.
Комментарии (17)
polybook
28.01.2016 17:34Какой-то очень суровый тест — неужели для решений диф.уравнений он используется?
Совпадение до 3-го знака (считал решение в середине стержня). Три расчета: для пространственной сетки в 100 и 200 точек (для 200 взял предыдущее и удвоенное число узлов по времени):
gbg
28.01.2016 21:16Это главный критерий получения зачета по численным методам, уж что-что, а точность должна соответствовать теории.
Увы, но теория разностных схем говорит, что знаков должно быть 4 (погрешность аппроксимации меньше квадрата шага сетки, который = 0.01). Так что вашей реализации НЕЗАЧТ.polybook
28.01.2016 21:21Буду признателен, если дадите ссылку на упоминание о 4-х знаках после запятой в теории разностных схем.
gbg
28.01.2016 22:43Вот например , раздел «Аппроксимация». Погрешность данной схемы составляет 2 порядка от шага по пространству (квадрат шага) и 1 порядок шага по времени: O(h^2+t). Из выполнения устойчивости по Куранту как раз следует, что t подобрана меньше h^2 в несколько раз, чтобы порядок аппроксимации по пространству и времени суммарно соответствовал.
Учебник Патанкара также могу рекомендовать.
Замечание в сторону — Технологию апостериорной проверки точности численных методов разработал Карл Рунге, зять того самого Куранта, который разработал условие устойчивости.polybook
29.01.2016 00:13Ваша первая ссылка ровно на явную схему бегущего счета, которая — один-в-один и рассмотрена в статье:
И о тесте про 4 знака после запятой там ничего нет.
Вообще-то моя статья — не про выбор хорошего алгоритма для решения диф.уравнения (явная схема Эйлера как раз плохая, а хороша она тем, что простая и часто работает). А про то, как быстро «на коленке» можно посчитать разностную схему в Маткаде, не требующем никаких навыков программирования (причем, в бесплатной версии).gbg
29.01.2016 00:30+2Дело не в конкретном числе «4 знака после запятой». Знаков схема дает не «4», а в зависимости от шага сетки. Если шаг сетки h<0.01, точность результата будет ± 0.0001 [O(h^2)]. (При условии, что длина = 1)
Обширный пласт теории численных методов вертится вокруг контроля точности вычислений. В противном случае, имеем классику компьютинга «на скриптах за 5 минут» — мусор на входе => мусор на выходе.
Такой результат не пройдет в качестве лабораторки в любом приличном учебном заведении — завернут. А по графикам вычислительные программы не тестируют.
Увы, но вычислительная математика крайне занудна. В качестве развлекухи для заманивания студентов она не подходит. Вот вы поставили 101 точку разбиения или 100?
Никто не будет искать, куда сбежала одна тысячная — график бабахнули и ладно. А потом спутники падают, мосты волнами идут, реакторы плавятся на АЭС, и так далее, и тому подобное.polybook
29.01.2016 00:46М=100 интервалов, это М+1=101 узлов.
Про занудность — в корне не соглашусь (моя миссия — как раз показать, что это не так). По мне, математика, как таковая (с формулами и доказательствами) — как раз скучна. А вот вычислительная математика — нет! Потому что я сразу чувствую связь с физикой и реальной жизнью. А, играя параметрами, могу проникнуть и в суть физических явлений, и в особенности вычислительных методов (как с устойчивостью схемы).
И конечно, в некоторых задачах допустимо считать «плюс-минус-лапоть», а в некоторых — нет. Для долговременных расчетов тех же координат ИСЗ (где обыкновенные диф.уры, но решения быстро осциллирующие, и точность критична) есть специальные алгоритмы.gbg
29.01.2016 01:05+1Вычислительная математика без доказательств никуда не годится, так как ее результаты ничем не обоснованы. Их можно применять в развлекательных целях, не более.
Вы реализуете вычислительный алгоритм, для которого есть гарантия по точности, вычисленная занудами у доски.
Вы эту гарантию не выполняете.
Напоминаю, что уменьшив шаг по пространству, нужно уменьшить шаг по времени, а количество по времени шагов — увеличить. Чтобы по времени полученные решения U(x,t) совпали.
Uranix
29.01.2016 09:47Ваше утверждение про точность в четвертом знаке неверно. Ошибка аппроксимации данной схемы —
, константа устойчивости — единица. Очевидно, что ошибка не фиксирована, она зависит от гладкости решения. Похоже, что о зануднсти вычислительной математики вы судите исключительно из своего опыта.
gbg
29.01.2016 10:05Моя оценка является более грубой, чем ваше E: h^2>E. То есть, на самом деле, схема чуточку точнее, и 4 знака давать она обязана.
Особенно, в такой постановке, когда есть аналитическое решение.
PapaBubaDiop
29.01.2016 00:36+1Ничо ни понял, для этой явной схемы курант надо брать 1/8 и все бу в шоколаде.
Astrei
28.01.2016 19:22Спасибо за статью, очень занимательно!
Вопрос немного не по теме. Раз уж речь зашла о тестах качества аппроксимации, не подскажите ли такой момент: для некой функции я провожу точечную аппроксимацию путем решения системы линейных уравнений. При этом я могу задать сколь угодно много уравнений и решить их, например, методом наименьших квадратов. Есть ли критерии, по которым я точно могу сказать, что достаточно трех/пяти/ста точек для успешной аппроксимации и рекомендации по выбору этих точек?polybook
28.01.2016 21:39Если посмотреть на самый простой случай — построение линейной регрессии Ах+В по массиву N точек х, у (т.е. одного уравнения), то легко получить оценку того, насколько близко расположены точки данных от линии регрессии. Будет зависимость от 1/(N-2)1/2. Подробнее см. в моей статье про корреляцию и регрессию (самый конец — последние три графика — как раз про такую оценку).
gbg
28.01.2016 22:47В общем виде таких рекомендаций нет и быть не может. Однако, тот же самый Рунге из комментария выше показал, что на равномерных сетках вас ждет Феномен Рунге, который заключатся в катастрофическом падении точности при увеличении количества узлов интерполяции полиномом.
gbg
Не хватает классических тестов качества аппроксимации — задать сетку из 101 и 201 точек, и при равном значении времени окончания эксперимента, убедиться в том, что точность полученного решения соответствует теории (должны быть верно вычислены (и должны совпасть) 4 знака после запятой, если метрическая длина стержня — 1 метр). Если этого нет — решение никуда не годится, увы.