Numerical simulation of physical fields by the collocation method.pdf Введение Одной из актуальных задач физики является численное моделирование полей различной физической природы. Наиболее часто для решения задач моделирования используются метод конечных элементов (МКЭ) и метод конечных разностей (МКР) [1, 2]. Однако при решении конкретных задач математической физики не всегда оправдано применение данных методов. Высокие требования к точности численного решения, необходимость учета граничных условий, содержащих производную по нормали от искомой функции на поверхностях со сложной конфигурацией, трудности, связанные с решением нелинейных задач, побуждают к поиску и разработке альтернативных численных методов решения краевых задач. Одним из перспективных численных методов решения задач математической физики является метод точечных источников поля (МТИ), который можно рассматривать как один из вариантов метода коллокаций [3-7]. В МТИ приближенное решение представляют в виде линейной комбинации базисных функций, точно удовлетворяющих основному уравнению. Однако в методе коллокаций допускается использование и других базисных функций, что значительно расширяет возможности численного метода [6-8]. Цель данной работы - развитие одного из вариантов метода коллокаций, который может быть использован при решении достаточно широкого круга краевых задач математической физики, в том числе и нелинейных задач. Метод полных коллокаций Пусть в области Ω двумерного или трехмерного пространства с границей ∂Ω решается краевая задача для уравнения , , (1) с условиями на границе ∂Ω , , (2) где , - линейные операторы в области Ω и на границе ∂Ω; и - заданные функции в области Ω и на границе ∂Ω. Пусть задана система линейно независимых функций , обладающих в области Ω свойством полноты. Представим решение задачи (1), (2) в виде ряда . Так как нахождение бесконечного числа коэффициентов разложения практически невозможно, то будем искать приближенное решение задачи (1), (2) в виде конечной суммы . (3) Для определения неизвестных коэффициентов разложения в (3) поступим следующим образом. Расположим в области Ω и на границе ∂Ω соответственно Nv и NG узловых точек, точек коллокации. Общее число точек коллокации N = Nv + NG должно соответствовать числу коэффициентов N в разложении (3). Координаты узлов, точек коллокации, в области Ω обозначим , а на границе ∂Ω . Для нахождения коэффициентов разложения ci подставляем (3) в уравнение (1) и в граничное условие (2) и потребуем их точного выполнения в точках коллокации. В результате получается система линейных алгебраических уравнений: ; (4) . (5) После решения системы линейных алгебраических уравнений (4), (5) получаем искомое численное решение в виде конечной суммы (3). Следует обратить внимание на то, что соотношение (3) дает приближенное аналитическое выражение для искомой функции . Это значит, что с полученным решением можно поступать как с любым другим аналитическим выражением. Его можно дифференцировать, производить другие действия, и при этом не возникает дополнительная численная погрешность. Описанный здесь вариант метода коллокаций, приводящий к системе линейных алгебраических уравнений в виде соотношений (4), (5), назовем методом полных коллокаций (МПК), а систему (4), (5) - системой МПК. Слово «полный» в названии МПК указывает на то, что требуется выполнение условий коллокации не только во всех внутренних узловых точках области решения Ω, но и на ее границе ∂Ω. Основная проблема, возникающая при нахождении приближенного решения (3), связана с возможностью плохой обусловленности системы (4), (5). Например, в [9] показано, что при интерполировании функции от двух переменных полиномами степени n определитель системы обращается в нуль, если узлы коллокации лежат на одной кривой порядка n. Аналогичный результат имеет место и для функций от трех переменных. Поэтому при регулярном, например, при равномерном, расположении точек коллокации избежать плохой обусловленности обычно не удается. Однако, если производить корректировку в расположении точек коллокации, производя случайным образом их небольшое смещение, то проблема, связанная с плохой обусловленностью системы МПК (4), (5), может быть устранена. Различные варианты методов коллокации отличаются друг от друга, в первую очередь, используемой системой базисных функций . Например, при решении двумерных задач в качестве базисных функций для МПК можно использовать функции вида , (6) а при решении трехмерных краевых задач использовать функции вида . (7) Каждому значению числа i в (6), (7) должен соответствовать определенный набор чисел ix, iy, iz. Кроме того, базисные функции необходимо упорядочивать таким образом, чтобы с увеличением порядкового номера базисной функции i суммарная степень базисной функции ni = ix+ iy+ iz не убывала, причем . Можно также искать приближенное решение в виде тригонометрической суммы (усеченного ряда Фурье). Определенную привлекательность имеют базисные функции, принимающие максимальное, например, единичное, значение в точке нахождения соответствующего узла коллокации и убывающие по мере удаления от этого узла. Это могут быть, например, функции следующих видов: ; (8) . (9) Здесь постоянные k подбираются таким образом, чтобы обеспечить наименьшую погрешность численного решения. Не исключается возможность использования других базисных функций для получения решения краевых задач с помощью МПК. МПК может успешно применяться для решения разнообразных задач математической физики. Ниже на конкретных примерах рассматриваются возможности МПК. Для этого приводятся результаты решения тестовых задач для эллиптических, параболических и волновых уравнений математической физики. Решаются задачи различной размерности: одномерные, двумерные и трехмерные, как линейные, так и нелинейные. Для каждого тестового примера решалась задача Дирихле и третья краевая задача. Наиболее естественным образом МПК применяется при решении краевых задач для линейных уравнений математической физики. Решение уравнений эллиптического типа с помощью МПК Ниже приведены результаты решения двумерной и трехмерной краевых задач. Ранее эти задачи решались при тестировании обобщенного метода точечных источников [10]. Сначала решалась двумерная краевая задача. В области Ω, представляющей собой квадрат со стороной L = 2, центр которого совмещен с началом координат, решалась задача для уравнения , , (10) с условиями на границе ∂Ω , . (11) Функции и в соотношениях (10), (11) подбираются таким образом, чтобы они соответствовали точному решению краевой задачи . Задавая целое число N0 и шаг по координатам X и Y в области Ω равным h = L/(N0-1), координаты точек коллокации вычисляются с помощью соотношений ; (12) , (13) где rnd - случайные числа из диапазона . Аналогичным образом задаются координаты узлов коллокации на границе ∂Ω, а также узлы коллокации при решении трехмерной краевой задачи. Приближенное решение краевой задачи (10), (11) искалось в виде полинома от двух переменных x, y, степень которого определялась полным числом точек коллокации N. В качестве базисных функций использовались степенные функции вида (6). Линейная независимость базисных функций (6) очевидна. При решении краевой задачи (10), (11) для каждого конкретного значения полного числа узлов коллокации N = N0 ∙ N0 с помощью соотношений (12), (13) определялось положение всех узлов коллокации, формировалась система линейных уравнений (4), (5), производилось численное решение задачи и оценивалась максимальная относительная погрешность решения в области Ω. Далее, при том же значении числа N решение задачи повторялось. Так как положение узлов коллокации, вычисленное с помощью соотношений (12), (13), при повторных решениях задачи случайным образом изменяется, то изменяется и значение максимальной относительной погрешности. В данной работе при каждом конкретном значении числа N краевая задача решалась 20 раз и определялись минимальная, максимальная и средняя из 20 значений максимальных относительных погрешностей. На рис. 1 представлена зависимость максимальной относительной погрешности от полного числа точек коллокации N, полученных при решении краевой задачи (10), (11). Кривая 1 на рис. 1 показывает зависимость средней (из 20) максимальной относительной погрешности от количества узлов коллокации, рассчитанная при использовании базисных функций (6). Видно, что с ростом числа N от начального значения N = 25 до значения N = 361 средняя максимальная погрешность решения уменьшается более чем на семь порядков, от εmax = 6.4∙10-2 до εmax = 2.7∙10-10. При дальнейшем увеличении числа точек коллокации наблюдается слабый рост погрешности до значения εmax = 9.5∙10-7 при N = 1089. Этот рост, видимо, связан с увеличением числа обусловленности системы линейных уравнений (4), (5), которое при N = 1089 имеет значение С = 7.0∙1019. Следует отметить, что диапазон вариации максимальной погрешности при каждом N не превышает двух порядков, т.е. в ряде случаев может считаться приемлемым. Рис. 1. Зависимость максимальной относительной погрешности от полного числа узлов N для двумерной задачи (кр. 1 и 2) и для трехмерной задачи (кр. 3) Краевая задача (10), (11) решалась также с использованием экспоненциальных базисных функций (8) (при k = 4). Кривая 2 на рис. 1 представляет зависимость погрешности этого решения от количества узлов N. Видно, что с ростом числа N от начального значения N = 25 до значения N = 721 погрешность решения уменьшается более чем на пять порядков, от εmax = 2.1∙10-1 до εmax = 1.4∙10-6 . Таким образом, при решении краевой задачи (10), (11) использование степенных базисных функций (6) дает более точное численное решение, чем при использовании экспоненциальных базисных функций (8). Разумеется, при решении других краевых задач более эффективным может быть использование экспоненциальных или иных базисных функций. Далее решалась трехмерная задача Дирихле. В области Ω, представляющей собой куб со стороной L = 2, центр которого совмещен с началом координат, решалась задача для уравнения , . (14) Функции в уравнении (14) и в правой части граничного условия подбирались таким образом, чтобы они соответствовали точному решению краевой задачи в виде . Кривая 3 на рис. 1 представляет зависимость относительной погрешности решаемой задачи Дирихле от размерности N, т.е. от полного числа уравнений в системе (4), (5). При вычислениях использовались степенные базисные функции (7). Видно, что уже при относительно небольшой размерности численной задачи N = 729 погрешность решения составляет εmax = 1.2∙10-8 . Решение нелинейных краевых задач с помощью МПК Моделирование нестационарных полей, удовлетворяющих уравнению параболического типа, может быть приведено к решению последовательности краевых задач для уравнений эллиптического типа. Для этого можно, например, аппроксимировать производную по времени от искомой функции разностным выражением вида , где τ - шаг по времени. Тогда, используя найденное решение , рассмотренным выше способом находится решение на последующем временном уровне . Если решаемая задача является нелинейной, то, при численном решении, нелинейности в уравнении или в граничных условиях следует вычислять на предшествующем временном уровне. Для решения нелинейных задач эллиптического типа используется метод установления [11]. В этом случае нелинейная краевая задача для эллиптического уравнения приводится к решению нелинейного параболического уравнения с последующей его линеаризацией и ростом времени для достижения стационарного режима. В качестве примера решалась тестовая задача Дирихле для нелинейного уравнения параболического типа (15) с известным точным решением [12] . (16) Так же как и в предыдущих примерах, решение находилось в области Ω, представляющей собой квадрат со стороной L = 2, центр которого совмещен с началом координат. Начальное и граничное условия для уравнения (15) определялись с помощью известного решения (16), которое использовалось также при вычислении погрешности численного решения. Решение краевой задачи и оценка его точности производились для момента времени t = 1. Определялась зависимость максимальной относительной погрешности от числа узлов коллокации для различных значений шага по времени τ. Результаты расчетов представлены на рис. 2. Рис. 2. Зависимость максимальной относительной погрешности от количества узлов при решении двумерной краевой задачи (15), (16) для нелинейного уравнения параболического типа На рис. 2 кривая 1 характеризует зависимость погрешности от количества узлов N при τ = 0.02, а кривая 2 - при τ = 0.08. Видно, что погрешность численного решения убывает с увеличением количества узлов коллокации. При этом погрешность при τ = 0.02 оказывается ниже, чем при τ = 0.08. Следует также обратить внимание на то, что при N > 200 дальнейшее увеличение количества узлов коллокации не ведет к заметному снижению погрешности. Это, видимо, связано с тем, что при N > 200 погрешность решения лимитируется шагом по времени, а не по координате. Из рис. 2 видно также, что при решении нелинейного уравнения параболического типа уже при N = 200 точность численного решения составляет доли процента, что может оказаться вполне приемлемым. Аналогичные результаты получаются при решении других, в том числе трехмерных, нелинейных задач. Ион молекулы водорода Важное значение имеют численные методы решения квантово-механических задач [8, 10, 13, 14]. МПК может успешно применяться при решении квантово-механических задач, задач на собственные значения. В качестве примера рассмотрим ион простейшей молекулы, молекулы водорода Н2+, состоящий из двух водородных ядер и одного электрона. Обозначим расстояние между двумя ядрами водорода (ядро А и ядро В) через а. Будем считать, что величина а изменяется адиабатически и при решении уравнения Шредингера ее можно принять постоянной. Это означает, что не будут учитываться колебания ядер около положения равновесия, а также вращение ядер вокруг центра тяжести. Движение электрона рассмотрим в полярной системе координат ( , z, ), поместив ядро А и ядро В на оси z на равном расстоянии от начала координат. Потенциальная энергия электрона в кулоновском поле ядер молекулы определяется выражением . Уравнение Шредингера для электрона с нулевым орбитальным моментом имеет вид , где Е - энергия электрона. Будем искать осесимметричные решения, для которых . Очевидно, такому решению соответствует, в частности, основное состояние иона. Воспользовавшись известным выражением для оператора Лапласа в цилиндрической системе координат, для волновой функции получим уравнение . Перейдем к атомным единицам длины и энергии . В этих единицах уравнение Шредингера имеет вид . (17) При численном решении задачи сначала полагали a = 2.1 ат.ед., что соответствует известному равновесному расстоянию между ядрами иона молекулы водорода [15]. Область решения Ω задавалась в виде квадрата со стороной L = 8a. Предполагалось, что на границах области Ω волновая функция обращается в нуль. Корректность такого предположения подтверждалась при последующем тестировании задачи. Система МПК (4), (5) для решаемой задачи является однородной и ее решение сводится к поиску таких значений энергии Е, для которых волновая функция имеет нетривиальное решение, обращающееся в нуль на границах области Ω. Для получения нетривиального решения задается некоторая точка М внутри области Ω, например, между ядрами водорода, в которой значение волновой функции заведомо не равняется нулю. В этой точке значение волновой функции полагается равным единице: . Это условие можно рассматривать как условие нормировки волновой функции. После получения решения волновую функцию можно легко перенормировать. Как показано выше, численное решение краевых задач с помощью МПК предполагает использование некоторой системы базисных функций. При решении уравнения Шредингера (17) обнаружилась невозможность использования степенных функций (6) из-за неустойчивости численного решения. Удовлетворительные результаты удается получить с помощью экспоненциальных (8) или гиперболических (9) базисных функций. Решение задачи производится в соответствии со следующим алгоритмом. Задаются начальное значение энергии электрона Е0 и шаг варьирования энергии δε. Желательно, чтобы этот шаг был много меньше расстояния ∆Е между ближайшими собственными значениями энергии, δε

Скачать электронную версию публикации