Размещено на http://www.allbest.ru/

Министерство образования и науки Российской Федерации

ФГБОУ ВО «Калмыцкий государственный университет имени Б.Б. Городовикова»

Магистерская программа «Математический анализ»

Факультет математики, физики и информационных технологий

Кафедра алгебры и анализа

ВЫПУСКНАЯ КВАЛИФИКАЦИОННАЯ РАБОТА

(магистерская диссертация)

«Численные методы в астрофизике»

Магистерская диссертация выполнена

студентом (кой) 2 курса

направления 01.04.01 Математика

Ольдеевой Екатериной Анатольевной

Научный руководитель:

кандидат физико-математических наук, доцент

Бисенгалиев Ренат Александрович/

Rdo.kz-2017



ОГЛАВЛЕНИЕ

ГЛАВА ПЕРВАЯ. МЕТОД НЬЮТОНА, МЕТОД РУНГЕ-КУТТА, МЕТОД СТРЕЛЬБ

1.1 Метод Ньютона

1.2. Метод Рунге-Кутта

1.3. Метод стрельб

ГЛАВА ВТОРАЯ. МОДЕЛЬ МАГНИТНОЙ АРКАДЫ НА СОЛНЦЕ. УЧЕТ ВЛИЯНИЯ ДИНАМИЧЕСКОГО ОХЛАЖДЕНИЯ ВЫСВЕЧИВАНИЕМ

2.1 Основные уравнения и равновесная модель

2.2 Постановка краевой задачи. Граничные условия

2.3 Алгоритм решения краевой задачи

2.4 Численный анализ закона дисперсии

ГЛАВА ТРЕТЬЯ. НЕУСТОИЧИВОСТЬ КОРОНАЛЬНЫХ ПЕТЕЛЬ

3.1 Модель корональной петли с продольным электрическим током

3.2 Линейные уравнения МГД в идеально проводящей среде

3.3 Обсуждение результатов



ГЛАВА ПЕРВАЯ. МЕТОД НЬЮТОНА, МЕТОД РУНГЕ-КУТТА, МЕТОД СТРЕЛЬБ

1.1 Метод Ньютона

Пусть оператор отображает в себя полное метричекое пространство X .

Рассмотрим итерационный процесс

(1.1)

решения уравнения

(1.2)

Определение. Если отображение удовлетворяет условию

(1.3)

при некотором и всех , то такое отображение называется сжатым.

Теорема. Если отображение - сжатое, то уравнение (1.2) имеет решение и

;

где .

Доказательство. Из уравнения (1.3), имеем



поэтому . При некотором имеем цепочку неравенств

(1.4)

Из критерия Коши имеем, что в последовательности есть некоторый предел X, переходя к которому при получаем

,

Очевидно, что

Отсюда , при n произвольном. В противном случае .

Теорема доказана.

Метод простой итерации будет иметь простую геометрическую интерпретацию тогда и только тогда, когда решается одно скалярное уравнение.

Компоненты приближений итерационного процесса в системе нелинейных уравнений можно найти из соотношений



(1.5)

…………………………………

Нахождение любого нового значения требует решения нелинейного уравнения с одним неизвестным

Переходный момент между итерационными методами (1.5) и (1.1) содержит метод, в котором компоненты приближений можно вычислить из соотношений

(1.6)

…………………………………

Пусть F(x) - оператор, который отображает линейное нормированное пространство X на линейное нормированное пространство Y, которое возможно совпадает с X. Если

(1.7)

при , где и -нормы в пространствах X и Y, то линейный оператор Р, действующий из пространства X в Y , называется производной оператора F(x).

Дальше будем обозначать оператор Р как F'(x). Пусть, например,

Если функции непрерывно дифференцируемы в окрестности заданной точки x, то

. (1.8)

Если за Р взять оператор умножения слева на матрицу

,

то совокупность соотношений (1.8) можно переписать в виде (1.7).

Оператор Р превращается в оператор умножения на производную в случае, если m=1.

В предположении существования производной , согласно (1.7), имеем

(1.9)

где X- решение уравнения , -некоторое приближение к Х.

При малой величине можно записать приближенное равенство

.

Так как F(x)=0, то



.

Пусть в приближение решения уравнения

оператор F' обратим , и его решение можно записать в виде

(1.10)

Данный итерационный процесс называется методом Ньютона.

Пусть .

Пусть при некоторых выполнены:

1. . (1.11)

2. (1.12)

при . Возьмем за

Теорема (о сходимости метода Ньютона). Итерационный процесс Ньютона (1.10) сходится при условиях (1.11) и (1.12) с оценкой погрешности

(1.13)

Доказательство. Пусть . С помощью индукции n докажем, что все . При некотором n и приближение , и если подставить в (1.12) получим следующее уравнение

.

Так как и F(x)=0, то соотношение может быть переписано в виде

Отсюда следует, что Из этого следует, что при для всех выполняется (1.14).

Пусть . Тогда после умножения на c неравенство (1.14) запишется в виде

Докажем справедливость неравенства с помощью индукции по n

.

При n=0 оно очевидно, поэтому предположим, что n=k, получим неравенство

Отсюда следует, что при всех n. И это означает, что

и отсюда следует (1.13). Из определения c и b, приближение при

.

Теорема доказана.

Метод Ньютона может измениться, если вычислить обращения оператора намного сложнее, чем значение . Изменения происходят следующим образом: выбирают или заранее задают некоторую возрастающую последовательность чисел .

Итерации производят по формуле

, при .

Увеличение числа итераций, сопровождающее часто такую модификацию, компенсируется большей «дешевизной » одного шага итерации. Выбор последовательности нужно производить с обоюдным учетом этих факторов.

В случае решения скалярного уравнения необходимо рассмотреть геометрическую интерпретацию метода Ньютона. Тогда расчетная формула (1.10) будет иметь вид

. (1.15)

При нахождении абсциссы точки пересечения с осью х касательной к кривой в точке , мы получим решение геометрически (рис.1).

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Рис.1 Рис.2

Также в случае, когда многочлен третьей степени, может произойти, что последовательность не сходится к корню при плохом начальном приближении. Например, на (рис.2) изображены все четные приближения совпадающие с а и нечетные - с b. В этом случае говорят, что метод «зациклился».

Если сравнивать асимптотическую скорость сходимости метода Ньютона и метода простой итерации, то согласно оценке погрешности

,

чтобы ошибка стала меньше ?, достаточно взять

.

В случае метода Ньютона правая часть (1.13) будет меньше ?, если

. (1.16)

Отсюда видим, что метод Ньютона, асимптотически, требует меньшего числа итерации при .

Также можно заметить, что в форме (1.10) метод Ньютона принимает вид разновидности метода простой итерации. Еще одну особенность метода Ньютона можно увидеть в случае скалярного уравнения , где производная правой части (1.15) по x равна . Отсюда, если , то и (рис.3) приобретёт вид (рис.4)



Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Рис.4



1.2 Метод Рунге-Кутта

Пусть на отрезке требуется найти решение дифференциального уравнения

(2.1)

при начальном значении , а аналитична в точке .

Проведя дифференциацию по х, получим следующие соотношения

Заменяя , а , получим значения

Отсюда, можно записать приближенное равенство

(2.2)

В случае, когда больше радиуса сходимости ряда

погрешность (2.2) не стремится к нулю при .

Иногда необходимо разбить на отрезки

.

Таким образом, по правилу: пусть значение уже найдено, вычисляем значения в точке производных решения исходного дифференциального уравнения, проходящего через точку , мы получим приближения к значениям решения

На отрезке полагаем

, (2.3)

то есть можем предположить, что

(2.4)

Рассмотрим метод Эйлера. Формула (2.3) при n=1 будет иметь вид

(2.5)

Для того чтобы, построить другой класс расчетных формул метода Эйлера необходимо: при известном значении решения вычислить значение .

Если рассматривать равенство

, (2.6)



то можно заметить, что если заменить интеграл правой части на величину , то погрешность будет иметь порядок , то есть

,

а так как , то будем иметь равенство вида

.

Если отбросить член порядка и заменить , , то мы получим уже рассмотренную расчетную формулу Эйлера (2.5). Эта формула не является точной, поэтому для того чтобы рассчитать точную формулу, мы должный аппроксимировать правую часть равенства (2.6) и использовать формулу трапеции, чтобы получить следующее равенство

.

В противном случае мы получим

(2.7)

Отсюда, точная расчетная формула Эйлера будет выглядеть следующим образом

(2.8)

Относительно , уравнение (2.8) будет неразрешимо. Поэтому будем дальше преобразовывать алгоритм.

В правой части равенства (2.7) заменим на некоторую величину

. (2.9)

Получим следующее равенство, в котором появляется новая величина , находящаяся между и :

При условии (2.9) величина имеет порядок ,а формула (2.7) примет вид

Из последнего соотношения можем определить следующие расчетные формулы

(2.10)

на том же шаге, с той же погрешностью и порядком, мы можем построить другую пару расчетных формул. Для этого нам необходимо заменить интеграл правой части равенства (2.6) по формуле прямоугольников

Или

Если , то

,

где есть результат вычисления по формуле Эйлера с шагом . Исходя из этих соотношений, получим пару расчетных формул

Эти методы являются частью семейства методов Рунге-Кутта. В процессе нахождения решений неизменны числа

Получаем следующий вид методов Рунге-Кутта



и можем предположить, что

Пусть . При достаточно гладкой функции своих аргументов , за гладкие функции параметра h будем принимать и . Сделаем следующее предположение при любых достаточно гладких функциях

а для некоторой гладкой

Отсюда, согласно формуле Тейлора, при где справедливо равенство

(2.11)

где - погрешностью метода на шаге, а s -порядок погрешности метода. При q=1 имеем

Далее за y будем принимать y(x). Из равенства выше можем сказать, что только в случае равенство выполняется при всех , а значению соответствует метод Эйлера. Для погрешности этого метода на шаге согласно (2.11), получаем выражение

В случае когда q=2, получим

где ,

Производные функции будут вычисляться в следующем виде:

,



Согласно исходному дифференциальному уравнению,

В выражения подставим значение h=0 и воспользуемся этими соотношениями; в итоге получится

(2.12)

Соотношение выполняется при всех f, если

, (2.13)

соотношение , если

и =0 (2.14)

Таким образом если выполнены три указанных выше соотношения ((2.13), (2.14)) относительно четырех параметров, то при всех f(x,y). Задавая один из параметров в любом порядке, получим различные методы Рунге-Кутта с погрешностью второго порядка малости по h. Для примера возьмем следующее, при получаем ,, что соответствует паре расчетных формул (10). При получем , что соответствует паре расчетных формул (11). Независимо от значений , , , , в случае уравнения y'=y, согласно (2.13), имеем . Тогда следует, что формул Рунге-Кутта со значениями q=2 и s=3 не существует.

Приведем результаты похожих вариантов для случая q=3. Расчетных формул, соответствующих значению s=4, не существует. Чтобы s=3, необходимо выполнить соотношения

, , ,

.

Эта система 6-ти уравнений с 8-ю неизвестными имеет бесконечное множество решений. Наиболее часто используемой совокупностью расчетных формул является

,

.

При q=4,5 не удается построить расчетных формул рассматриваемого вида со значением s=5.

При q=4 имеется двухпараметрическое множество расчетных формул, соответствующих s=4. Для примера возьмем одно параметрическое семейство подобных формул:

,

.



Наиболее часто используемой совокупностью формул этого семейства можно считать ту, которая соответствует t=1:

(2.15)

Среди других совокупностей формул со значениями q=s=4 отметим совокупность

(2.16)

При дальнейшем рассмотрении будет важно, что погрешность метода на шаге имеет главный член, а именно, справедливо представление вида

. (2.17)

Представим основные этапы доказательства этого соотношения. Предположим, что правая часть и все ее производные до порядка s+1 включительно ограничены равномерно в области G:

.

Отсюда также будут равномерно ограничены производные всех решений уравнения до порядка s+2, включая и его самого. Согласно формуле Тейлора, равенство (2.11) можно записать в уточненной форме

Будем иметь равенство

Нетрудно заметить, что величины явно выражаются через значения в точке (x,y) функции f и ее производных порядка не выше s. Примеры таких явных выражений мы уже получали.

Поскольку правая часть дифференцируема s+1 раз, то отсюда следует, что функция дифференцируема в области G и ее производные равномерно ограниченны в этой области. Аналогично устанавливается, что величина равномерно ограничена при . Некоторые частные случаи этого утверждения также следуют из выписанных выше соотношений.

1.3 Метод стрельб

Численный метод, заключающийся в сведении краевой задачи к некоторой задаче Коши для той же системы дифференциальных уравнений, называется метод стрельбы. Для того чтобы рассмотреть этот метод, возьмем пример простой задачи для системы двух уравнений 1-го порядка с краевыми условиями достаточно общего вида

, , , (3.1)

, . (3.2)

Для начала выберем любое значение , и рассмотрим левое краевое условие каждого в виде алгебраического уравнения . Далее определим значение , соответствующее этому уравнению. Используем значения , в качестве начальных условий задачи Коши для системы (3.1). Затем эту задачу нужно провести через интеграцию любым численным методом. В результате получим решения , , зависящих от , как от параметра.

Найденное решение будет удовлетворять левому краевому условию (3.2) при правильно подобранном значении . И все же правому краевому условию это решение, вообще говоря, не удовлетворяет: при его подстановке левая часть правого краевого условия, рассматриваемая как произвольная функция параметра :

, (3.2)

не будет равно нулю. Требуется любым возможным способом изменять параметр , до тех пор, пока не подберем нужное значение, для которого с заданной точностью. Отсюда следует, решение краевой задачи (3.1) свелось к нахождению корня одного алгебраического уравнения

. (3.3)

Рассмотрим, различные методы решения этой задачи, которые удобнее всех применять в данном случае.

Наиболее простым является метод дихотомии. В нем производят пробные « выстрелы » -- расчеты с наудачу выбранными значениями , до того момента, как среди величин не окажется различных по знаку. Пара таких значений , образует «вилку». Последовательно деля ее пополам до нужной точности, производим «пристрелку» параметра . С помощью этого процесса весь метод получил название стрельбы.

Однако нахождение каждого нового значения функции достаточно трудоемко, т.к. требуется численно интегрировать системы (3.1a). В таком случае корень уравнения (3.3) желательно находить более быстрым методом, чем дихотомия.

будет иметь непрерывную производную, в том случае, если правые части уравнений (3.1a) и левые части краевых условий (3.1б) имеют непрерывные и ограниченные первые производные. В итоге можно построить аналог метода Ньютона. Пока нам известен только способ вычисления и нужно также научиться определять производную

. (3.4)

Производные по параметру от решения задачи Коши входящие в (3.4) можно найти, если продифференцировать по этому параметру систему (3.1a).

Обозначая

, (3.5)

и дифференцируя (3.1a) по параметру, получим

, (3.6a)

,

Одно из начальных условий для этой системы очевидно: ; второе условие нетрудно найти, продифференцировав левое краевое условие (3.1б) по . Отсюда получаем

, . (3.6б)

Для того чтобы, определить вспомогательные функции , , проинтегрируем систему (3.6а) с начальными условиями (3.6б) совместно с задачей Коши для системы (3.1а). И подставив их значения при в (3.4), найдется значение производной правого краевого условия по пристрелочному параметру. По формуле касательных определим новое значение параметра:

. (3.7)

Однако описанный выше способ требует интегрирования лишней пары дифференциальных уравнений, что приводит к усложнению и двукратному увеличению трудоемкости каждой итерации. Поэтому им пользуются не часто.

Если решать уравнение (3.3) разностным аналогом метода Ньютона т. е. методом секущих, можно избежать этого усложнения. Для этого первые два расчета делают с наудачу выбранными значениями , , а последующие значения параметра вычисляют по формуле:

. (3.8)

Вместо этого процесса можно использовать метод парабол, в котором также не требуется располагать явным выражение производных, а достаточно лишь знать об их существовании. Напомним, что последние три метода быстро сходятся вблизи корня; сходимость вдали от корня зависит от того, насколько удачно выбрано нулевое приближение.

Линейные задачи решаются методом стрельбы особенно просто. Пусть система (3.1а) и краевые условия (3.1б) линейны;

, , (3.9а) ,

. (3.9б)

Тогда начальные условия соответствующей задачи Коши примут вид

, . (3.9в)

Легко понять, что решение задачи Коши (3.9a), (3.9в) будет линейно зависеть от параметра , поэтому также будет линейной функцией. Но линейная функция одного аргумента полностью определяется своими значениями в произвольных двух точках и , а ее график является прямой, т. е. совпадает со своей секущей. Значит, найденное по формуле секущих (3.8) значение является точным корнем уравнения (3.3), так что расчет с этим значением параметра даст искомое решение. Таким образом, для решения линейной краевой задачи (3.9а)--(3.9б) достаточно трижды решить задачу Коши.

Замечание. Для линейных задач можно немного уменьшить объем расчетов, если воспользоваться тем, что общее решение линейной неоднородной системы равно сумме ее какого-нибудь частного решения и общего решения соответствующей однородной системы. Найдем частное решение неоднородной системы (3.9а), (3.9в), которое соответствует значению , и обозначим его через , . Далее рассмотрим соответствующую однородную задачу Коши

, ,

, ;

вычислим ее решение и обозначим его через , . В этом случае общее решение неоднородной задачи Коши, удовлетворяющее левому краевому условию (3.9б) является однопараметрическим семейством

, . (3.10)

Значение параметра выбираем так, чтобы удовлетворить правому краевому условию (3.9б):

.

Затем, чтобы избежать третьего интегрирования задачи Коши, найдем искомое решение по формуле (3.10).

Метод стрельбы прост, т. к. применим как к линейным, так и к нелинейным задачам и позволяет использовать при численном интегрировании схемы Рунге -- Кутта высокого порядка точности. К большинству задач типа (3.1) он применяется успешно.

Если краевая задача (3.1) хорошо обусловлена, а соответствующая ей задача Коши плохо обусловлена, возникают некоторые затруднения. При этом численное интегрирование задачи Коши определяет функцию с большой погрешностью, что осложняет организацию итераций.

Таким образом, для данного случая пробуют поставить начальные условия на другом конце отрезка , т. е. интегрировать задачу Коши справа налево; в большинстве случаях устойчивость улучшается. Если изменение направления интегрирования не помогает, то такую краевую задачу решают, либо специальными, либо разностными методами.

Одним из специальных методов для линейных краевых задач является дифференциальная прогонка. Этот метод хорошо устойчив именно в том случае, когда задача Коши для исходной линейной системы плохо обусловлена. Однако при хорошей устойчивости линейной задачи Коши прогонка становится недостаточно устойчивой. Поэтому в настоящее время дифференциальная прогонка употребляется не часто.

При исследовании нестационарностей и, в частности, волновых движений в газовых и плазменных подсистемах астрофизических объектов часто приходится применять линейный анализ устойчивости.

При решении задач об устойчивости прежде всего определяются равновесные стационарные распределения , являющиеся решением системы уравнений гидродинамики или магнитной гидродинамики (МГД), в которых полагается (здесь радиус вектор точки рассмотрения, время, - любой из термодинамических параметров, компоненты скорости среды или магнитного поля в МГД - случае).

Далее применяется стандартная процедура линеаризации, для чего все физические величины, характеризующие состояние рассматриваемого объекта, представляются в виде: , где , и производится пренебрежение квадратичными по малым возмущениям слагаемыми (с математической точки зрения это отвечает разложению в ряд Тейлора всех входящих в систему уравнений гидродинамики функции с последующим учетом лишь линейных по малым приращениям слагаемых; при этом слагаемые, не содержащие таких приращений, очевидно, взаимно сокращаются в силу уравнений стационарного баланса).

Следующим шагом, как правило, является применение к полученной линейной системе дифференциальных уравнений в частных производных метода нормальных мод с целью сведения ее к системе двух обыкновенных дифференциальных уравнений.

Метод нормальных мод заключается в следующем.

Хорошо известно, что произвольное возмущение термодинамических параметров и компонент скоростей можно разложить в ряд Фурье по гармоникам вида:

. (3.11)

При этом и связаны законом дисперсии (зависимость частоты от волнового вектора). В случае, если закон дисперсии допускает не дискретный, а непрерывный спектр частот, сумма заменяется интегралом.

Однако если речь идет о неустойчивых возмущениях, то экспоненциально нарастающая во времени гармоника с максимальным инкрементом (инкремент ) станет доминирующей для определения структуры течения, и об остальных гармониках можно просто забыть, так как они не дадут существенного вклада в морфологию рассматриваемой системы.

Поэтому при решении задач о линейной устойчивости гидродинамических систем выявляются исключительно гармоники, обладающие максимальным инкрементом при данных значениях параметров системы, а остальные гармоники игнорируются. В связи с этим индекс «» мы далее не используем.

В принципе это напоминает хорошо известный и широко используемый в радиофизике метод комплексных амплитуд, но, естественно, при обсуждении физики неустойчивостей необходимо выбирать действительную часть от комплексной амплитуды возмущений.

Довольно часто симметрия задачи оказывается такой, что коэффициенты линеаризованной системы уравнений гидродинамики в частных производных зависят только от одной пространственной координаты; для определенности полагаем, что от декартовой системы координат. Тогда, в силу однородности (независимости) коэффициентов системы по координатам , и , решение ищется в виде:

. (3.12)

Подстановка решения вида (3.12) позволяет разделить переменные и свести систему уравнений в частных производных к системе двух обыкновенных дифференциальных уравнений на описывающие комплексные амплитуды возмущений функции и . Как правило, в качестве таких функций выбираются амплитуды возмущенного давления и возмущенного лагранжева смещения в направлении такого, что , где - амплитуда возмущенной компоненты скорости. В наиболее общем случае такая система выглядит следующим образом:

, (3.13)

. (3.14)

Здесь , , и - комплексные коэффициенты, зависящие от стационарных (равновесных) распределений термодинамических параметров модели, скорости и магнитного поля, от компонент волнового вектора , и комплексной частоты . Вместе с граничными условиями на собственные функции и , определяемыми из физических соображений, система (3.13)-(3.14) образует краевую задачу типа Штурма Лиувилля на собственные значения частоты .

Как правило, для имеющих практический научный интерес задач аналитические решения получить невозможно - это уже сделали задолго до нас. Чаще всего поставленную краевую задачу приходится решать численно на ЭВМ методом стрельб; однако, тем не менее, в ряде наиболее простых случаев она имеет аналитические решения и сводится к дисперсионному уравнению (степенному алгебраическому или трансцендентному). Именно такие простые случаи используются для тестирования разрабатываемых компьютерных программ, о чем пойдет речь позднее.

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

, , , , .

Здесь - расстояние между границами, - адиабатическая скорость звука в среде, - равновесная плотность среды, - стационарная скорость среды.

Естественно, при операции обезразмеривания необходимо соблюдать правила математики, чтобы не нарушить тождеств. Например, чтобы обезразмерить левую часть уравнения (3.13), поступаем следующим образом:

.

Мы однако, чтобы не загромождать изложение, сохраним прежние обозначения, подразумевая, что все переменные уже обезразмерены.

Метод стрельб получил свое название из-за хорошо просматривающейся аналогии со стрельбой из артиллерийского орудия; перелет-недолет, отклонения вправо-влево. Задача наводчика найти такие углы возвышения и горизонтальной наводки, при которых снаряд попадает в цель, причем сделать это максимально быстро.

Применительно к обсуждаемой системе уравнений (3.13)-(3.14), для использования метода стрельб необходимо проделать операции, описанные ниже.

1. Из физических соображений определить значение функции и на границах и рассматриваемого пространственного интервала. Причем сделать это нужно, по возможности, точно и корректности полученных результатов. Аналогом значений , , и являются местоположения артиллерийского орудия и цели соответственно.

2. Задать наиболее правдоподобное, с Вашей точки зрения, начальное значение комплексной частоты , действительная и мнимая части которой будут соответствовать, образно говоря, углам возвышения и горизонтального наведения (чтобы, по крайней мере «стрелять» в сторону цели, а не в противоположную).

3. Произвести «выстрел», т.е. проинтегрировать систему (3.13)-(3.14) от до . Для обсуждаемого типа задач это обычно делается с использованием метода Рунге-Кутта четвертого порядка точности. Распределения и символизируют при этом траекторию полета снаряда.

4. Продолжать «огонь», изменяя «прицел», т.е. значения частоты таким образом, чтобы с -го выстрела «снаряд» попал в цель (выполнились граничные условия при ). В нашем случае роль «наводчика» будет исполнять итерационный метод Ньютона, а будет номером итерации, при котором выполнились указанные граничные условия.

Идеология и алгоритм метода Рунге-Кутта подробно описаны выше. Заметим лишь, что метод Рунге-Кутта так же корректно работает с комплекснозначными функциями, как и с действительными - нужно лишь описать в программе эти функции и комплексные переменные типом COMPLEX*16 (число обозначает использование двойной точности для работы с переменными).

Иначе обстоит дело с методом Ньютона. Хотя описание этого метода есть во всех учебниках по численным методам, и его алгоритм формально выглядит предельно просто, при работе в комплексной области возникают «подводные камни», с которыми авторы неоднократно сталкивались. Вот на этом и хотелось бы остановиться поподробнее, чтобы сэкономить читателям время, силы и нервы.

Формула алгоритма метода Ньютона действительно проста; пусть требуется найти численное решение уравнения , тогда:

, (3.15)

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

В областях значений параметров задачи, для которых метод работает стабильно, модуль последнего слагаемого в (3.15) быстро уменьшается, т.е. сходится с заданной точностью к истинному значению корня за 3-5 итераций. Это очень быстрая сходимость по сравнению с другими итерационными методами (дихотомии, простой итерации и т.д.) для которых нормой являются 15-30 итераций. Кроме того, алгоритмы указанных методов принципиально неприменимы для работы с комплекснозначными функциями.

Однако, к сожалению, метод Ньютона не только не обладает абсолютной сходимостью, но, и наоборот, очень чувствителен к поведению функции вблизи истинного корня и к близости “угаданного” начального приближения к этому корню. В определенных ситуациях вместо того, чтобы последовательность стягивалась к истинному решению с ростом , она уходит от него все дальше и дальше! В таких случаях говорят, что метод расходится. Условие сходимости метода Ньютона можно найти, например, в книгах Калиткина и Корн, Корн (1974), однако, поскольку мы работаем с комплексными числами, оно несколько видоизменяется, относительно приведенного в книге Калиткина, и выглядит так:

<. (3.16)

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

1. Используйте вычисления только не двойной точности - накопление погрешности только ухудшит сходимость метода Ньютона.

2. Не задавайте слишком высокую точность в условии прекращения итераций, когда можно считать, что корень успешно найден. Вполне достаточно

3. Поставьте ограничитель - если за тридцать итераций корень не найден, смело выходите из цикла и прерывайте процесс, иначе компьютер «уйдет в себя» и будет бесполезно трудиться, пока не возникнет переполнение.

4. Вычисляете численно - при аналитических вычислениях велик шанс сделать ошибку в математических выкладках, тем более, что в реальных задачах формулы достаточно громоздки. Кроме того, в силу последней причины может возникнуть накопление ошибок, если в аналитическом выражении для производной придется отнимать или складывать большие числа и малые. Тогда, желая повысить точность работы метода Ньютона. Вы наоборот, ее понизите. Вполне достаточно взять производную по двум - трем точкам.

5. Причисленной аппроксимации производной не повторяйте распространенной ошибки вида:

(3.17)

Это очень сильно ухудшает сходимость метода даже в «здоровом» диапазоне параметров, где при правильном дифференцировании он бы прекрасно работал. Задайте в самом начале программы один шаг дифференцирования и пользуйтесь только им, применяя формулу:

(3.18)

6. Не выбирайте модуль шага слишком маленьким, это также сильно ухудшает сходимость - если поверхности действительной и мнимой частей функции будут достаточно пологими относительно плоскости, в числителе выражения (3.18) будет разность очень близких чисел, которые, не будем забывать, вычисляются с погрешностью. Может случиться ситуация, когда действительно значащие цифры взаимно сократятся и вместо производной мы получим абсолютно случайное число. Вполне достаточно, чтобы выполнялось конкретный выбор шага зависит от обусловленности задачи. Часто студенты (и не только студенты) спрашивают, в каком направлении в плоскости лучше выбирать это приращение вдоль действительной оси, вдоль мнимой, или под равным углом к обеим? В действительности это совершенно безразлично и не влияет на сходимость метода Ньютона (по - меньшей мере, один из составителей пособия это тщательно проверял и потратил на это немало времени).

7. Часто в литературе встречается утверждения, что если выполняется условие: , то итерации можно прекращать и выходить из цикла. Это неверно. То, что мы подошли близко к точке корня, еще не означает, что , а ведь именно этого требует исходное уравнение. И наоборот. Это легко понять на примере действительной функции одного переменного - поскольку метод Ньютона является методом касательных (за значение принимается точка пересечения с осью абсцисс касательной к графику функции, проведенной в точке ), - то для очень пологих функций, пересекающих ось абсцисс под малым углом, значение модуля функции может быть очень малым, но разброс ее аргументов, вычисленных на последующих итерациях - большим. И наоборот, для крутых функций, графики которых пересекают указанную ось под углом, близким к вертикали, разброс указанных аргументов может быть малым, а значения модулей функции для этих значений аргументов - большими. Читателям предлагается вооружиться бумагой и ручкой и убедиться в это лично, проитерировав по методу касательных вручную на графике.

Поэтому мы, основываясь на личном опыте, предлагаем использовать смешанное условие для обсуждаемой нами краевой задачи типа Штурма - Лиувилля, решаемой нами методом стрельб: если или , то возврат в цикл для продолжения итераций (здесь индексами «» помечены граничные значений функций, которые по аналитическим соображениям должны иметь в точке ).



ГЛАВА ВТОРАЯ. МОДЕЛЬ МАГНИТНОЙ АРКАДЫ НА СОЛНЦЕ. УЧЕТ ВЛИЯНИЯ ДИНАМИЧЕСКОГО ОХЛАЖДЕНИЯ ВЫСВЕЧИВАНИЕМ

Развитие волноводно-резонансной центробежной МГД - неустойчивости во вращающемся цилиндрическом слое замагниченной плазмы может привести к изменению этого слоя к виду, похожему на солнечные магнитные аркады, при этом не беря в учет эффект динамического охлаждения высвечиванием. Тем временем, можно увидеть наличие сильного высвечивания на стадии разделения рассматриваемого слоя на отдельные арки. Далее покажем, что учет охлаждения высвечиванием усиливает неустойчивость вращающегося цилиндрического слоя замагниченной плазмы, ускоряя его разделение на аркадную структуру. Будем работать в границах линейного анализа устойчивости. Моделью волокна с вращательными движениями плазмы будем считать цилиндрический слой, скорость и магнитное поле которого имеют лишь азимутальную часть. При этом в остальном пространстве магнитное поле устремлено вдоль образующей цилиндра, а плазма будет находиться в состоянии покоя.

Потому как в обычной ситуации нижняя часть этого цилиндра опущена в фотосферу, где он не иначе как перестанет быть цилиндром, наше рассмотрение применимо или для коротковолновых возмущений в азимутальном направлении (при больших номерах индексов ), или для осесимметричных возмущений, не имеющих азимутальной структуры ().

Поскольку мы рассматриваем устойчивость значительно упрощенной стационарной модели, будем также использовать интегральную функцию охлаждения , где Т- абсолютная температура среды.



2.1 Основные уравнения и равновесная модель

Исходной системой уравнения МГД считают

(1.1)

(1.2)

(1.3)

(1.4)

и уравнения состояния

(1.5)

где Р - давление, г - коэффициент адиабаты, е - плотность внутренней энергии, R- универсальная газовая постоянная, м - молярная масса вещества, Т - абсолютная температура среды, - квадрат адиабатической скорости звука, S - удельная энтропия единицы объема среды.

Из уравнения (1.1) мы принимаем во внимание вероятность охлаждения газа высвечиванием (q<0) или его нагрев внешним излучением (q>0), где - количество энергии, используемое единицей массы вещества в единицу времени, - функция нагрева, - функция охлаждения, зависящие лишь от температуры Т. Эффект теплопроводности игнорируем, так как он намного меньше двух других указанных.

Из уравнения (1.1) вычтем (1.4), которое в свое время умножено на , также вычтем из (1.1) уравнение Эйлера (1.2), домноженное на сV, и уравнение индукции (1.3), умноженное на В/4р, далее, с помощью тождества векторного анализа, получим уравнение только на внутреннюю энергию:

(1.6)

Вычислим уравнение для замыкания системы (1.2)-(1.6), используя при этом уравнение состояния и полную производную по времени:



(1.7)

где - уравнение состояния идеального газа, и переменны с индексом «0» являются начальными равновесными стационарными значениями величин.

То отсюда, последнее уравнение системы принимает вид:

(1.8)

Примем во внимание, что здесь и далее , где r - радиус в цилиндрической системе координат.

В качестве модели будем использовать цилиндрический слой, размещенный горизонтально над поверхностью верхней фотосферы (рис.2.1), вращающаяся с постоянной угловой скоростью Щ.



В слой между внутренним радиусом и наружным магнитное поле и скорость вмороженное в это поле вещества V азимутальны, а в остальном пространстве магнитное поле направленно вдоль образующей цилиндрического слоя, и плазма изначально не подвижна. То отсюда, начальное равновесное распределение параметров принимает вид:

(1.9)

Равновесную плотность вещества внутри и снаружи считаем постоянными, где индексом «in» обозначается параметры при , а индексом «ех» параметры при и .

Для этой модели можно игнорировать воздействие гравитации, так как оценки указывают на то, что центробежная сила может быть на много больше. корональный петля магнитный гидродинамика

На тангенциальных разрывах должно выполняться классическое условие равенства полных давлений:

(1.10)

при .

Стационарность модели определяется условием равенства линейной скорости вращения и альфвеновской скорости, то есть



.

Рассмотрим МГД-уравнение Эйлера:

(1.11)

Рассматривая стационарное состояние, получим:

(1.12)

или

(1.13)

Используя (2.10), находим:

(1.14)

Решая полученное равенство в цилиндрической системе координат, учитывая, что , получим:

(1.15)



Также в данной модели выполняются равенства:

(1.16)

(1.17)

Для осуществления стационарного баланса должно выполняться условие . Так как при выполнении линейного анализа устойчивости перемена температуры проходит в узком диапазоне вблизи его равновесного значения, почти всегда зависимость для межзвездной среды с хорошей степенью точности можно аппроксимировать степенными функциями. Предположим, что выполняется , где - нормировочная постоянная. Аналогичным образом полагаем, что .

То есть это значит, что выполняется

. (1.18)

Применяя стандартную процедуру линеаризации, мы представляем все величины в виде суммы равновесного значения и малого возмущения :



(1.19)

при ,

при .

Далее выпишем линеаризованные системы уравнений МГД для цилиндрического слоя(), и он будет иметь вид:

Решение ищем в виде

Отметим, что выполняются равенства

Подставим эти равенства в систему:



Подставляем вместо



Используя силу однородности коэффициентов линеаризованных уравнений по ц, z и t, мы можем найти решение в виде где -комплексная амплитуда, -комплексная частота возмущений.

В таком случае, линеаризованная система сводится к системе из двух о.д.у вида:

(1.21)

(1.22)

где -возмущенное радиальное лагранжево смещение среды, которое связанно с радиальной компонентой скорости соотношением:

(1.23)

где - частота с учетом доплеровского сдвига,



Используя эти обозначения, перепишем стационарный баланс в безразмерном виде:

(1.24)

где - альфвеновская скорость в слое, - альфвеновская скорость во внешней среде (относительно слоя).

Из уравнения (1.13) следует:

(1.25)

где М-число Маха.

Вне цилиндрического слоя ()линеаризованная система МГД принимает вид:

(1.26)

Используя силу однородности коэффициентов линеаризованных уравнений по ц, z и t, мы можем найти решение в виде

Введем обозначения:

(1.27)

Система (1.27)-(1.28) сводится к видоизмененному уравнению Бесселя:

(1.29)



где

(1.30)

- комплексная амплитуда возмущенного лагранжева смещения во внешней среде.

2.2. Постановка краевой задачи. Граничные условия

Используя условия необращения амплитуды возмущений в бесконечность при , получим:

(2.1)

(2.2)

где -амплитудные постоянные, - видоизмененная функция Бесселя первого рода, - функция Макдональда.

Помимо этого, необходимо выполнение непрерывности возмущенного ридиального лагранжева смещения на изогнутых границах тангенциальных разрывов при :

(2.3)

Существование динамического охлаждения высвечиванием меняет баланс возмущенных сил при :

(2.4)

где

(2.5)

(2.6)

То отсюда краевая задача будет иметь вид:

1)система уравнений в слое:

(2.7)

(2.8)

2) вне слоя:



(2.9)

(2.10)

2.3 Алгоритм решения краевой задачи

Поставленная краевая задача решалась численно.

Алгоритм решения краевой задачи:

1) Задаем начальное значение , где - безразмерная ная частота;

2) Из (2.1) вычисляем . Предположим, что , так как при рассмотрении линейного анализа устойчивости берем любой выбор безразмерных ных постоянных;

3) Определим из уравнения . Для этого изпользуем известную формулу для первой производной модифицированной функции Бесселя 1-го рода:

. (2.11)

Используя формулу (2.3), мы можем определить

4) Так как вы уже вычислили значения ,, , поэтому, в силу (2.4), мы можем также вычислить ;

5) Интегрируем систему уравнений (1.21) - (1.22) с помощью метода Рунге-Кутта от до ;

6) Определим , используя (2.3);

7) Найдем из уравнения . Для этого используем формулу для первой производной модифицированной функции Бесселя 2-го рода:

, (2.12)

и из (2.3) можем определить ;

8) Стремимся к выполнению граничного условия (2.4) при . Для того чтобы (2.4) выполнилось подставим значения ,, , и применим итерационный метод Ньютона.

2.4 Численный анализ закона дисперсии

Ранее было сказано, что поставленная краевая задача типа Штурма-Лиувилля решалась численно на ЭВМ методом стрельб. Основные результаты численных расчетов были получены при следующих значениях термодинамических параметров вещества: плотности и , магнитных полей и , скорости звука и , температуры и . Для значений температур показатель степени охлаждения изменяется в основном от -1 до 0 при . Для функции нагрева чаще всего выполняетя , где С - удельная теплоемкость газа, - его начальная и текущая температуры, поэтому .

Из рис.2.2 видно для осесимметричных возмущений (), развиваются 4 семейства неустойчивых мод (рис. 2.3). Гармоники каждого семейства имеют ярко выраженный волноводно-резонансный характер, поскольку различны друг от друга числом нулей возмущенного лагранжева смещения между границами слоя. Для интерпритации их физической природы полезно привлечь асимптотические решения.

Рис. 2.2. Зависимости безразмерных частоты(слева) и инкремента осесимметричных (m=0) неустойчивых возмущений (справа) в единица угловой скорости вращения от безразмерного волнового числа. Цифрами обозначено число нулей безразмерного возмущенного лагранжева смещения по радиусу в слое. M=1.5, s=6.0, F=0.475, 1.0, c=-0.5, h=1, Штрих-пунктирными линиями показаны асимптотики: I-быстрых магнитозвуковых волн в слое (), II- медленных магнитозвуковых волн во внешней относительно слоя среде.

Рис. 2.3. То же, что на рис.2.2, для ММЗВ - энтропийных волн

Локальный закон дисперсии допускает четыре типа решений для возмущений, распространяющихся вдоль поля во внешней относительно слоя среде:

(2.13)

знак «+»под радикалом отвечает быстрым магнитозвуковым волнам (БМЗВ), «-» - медленным магнитозвуковым волнам (ММЗВ), распространяющимся вдоль аркады в z-направлении.

Для фазовой скорости БМЗВ и ММЗВ в слое будет отсутствовать последнее слагаемое внутри слоя для оссесиметричных ( m=0) возмущений в выражении потому, что в него входит компонента поля, направленная вдоль волнового вектора, то есть получим

(2.14)

Аналогично найдем для частоты

(2.15)

Тогда получим

(2.16)

Таким образом, для БМЗВ имеем

(2.17)

а для ММЗВ частота оказывается вырожденной и равной частоте энтропийных возмущений внешней среды для адиабатического () случая:

(2.18)

Также можно предположить о наличии гироскопических мод вращающегося слоя с частотами:

(2.19)

Для первого семейства (рис.2.2) мы идентифицируем их как гироскопические моды так, как основная безусловная мода сходится к Они обусловлены вращением плазмы в аркадах. Также известно, что для гироскопических мод присуще слабая зависимость частоты от в длинноволновой области.

Для второго семейства, мы идентифицируем гармоники как возбуждаемые ММЗВ, которые распространяются вдоль слоя снаружи от него при выполнении асимптотики ММЗВ (2.13).

Для третьего семейства мод, при правильной работе асимптотики (2.17), мы разделяем их на классы как гармоники БМЗВ, распространяющихся внутри слоя. Также аргументом является уменьшение действительной части их частоты с ростом внешнего поля, чего не наблюдалось бы для БМЗВ внешней среды (2.13).

С определение четвёртого генезиса семейства волноводных гармоник (дисперсионные кривые, которого приведены на рис.2.3) обстоят дела намного сложнее. По виду их частота близка к нулю, то есть близка к частоте ММЗВ внутри слоя, а также близка и к частоте энтропийных возмущений внешней среды. Мы можем сделать вывод, что гармоники четвертого семейства имеют смешанный, гибридный характер ММЗВ - энтропийных волн, поскольку в одном случае, инкремент не обращается в ноль в отсутствии высвечивания, хотя оно должно быть для энтропийных волн, а в другом случае этот инкремент сильно нарастает с увеличением интенсивности высвечивания, что как раз и свойственно для энтропийной моды.

С учетом условия (1.15) для фазовой скорости БМЗВ и ММЗВ в цилиндрическом слое имеем

(2.20)

Для частоты возмущений, которые распространяются вдоль поля, с учетом доплеровского сдвига получаем

(2.21)

Тогда

(2.22)

Можем сказать, что для БМЗВ выполняется

(2.23)

а для ММЗВ

(2.24)

Если , то для того чтобы получить частоту, необходимо умножить на . В итоге получим

(2.25)

где в качестве используется значение .

Соответственно, для БМЗВ имеем

(2.26)

а для ММЗВ имеем

(2.27)

Асимптотические решения (2.23)-(2.24) и (2.19), которое отвечает гироскопическим модам, показаны на рис. 2.4 штрих-линиями.



Рис.2.4. Зависимости безразмерных частоты (слева) и инкремента неустойчивых возмущений (справа) в единицах угловой скорости вращения для азимутальной моды m=6 от безразмерного волнового числа. Остальные параметры такие же как на рис. 2.2. Штрих-линиями показаны асимптотики: I - асимптотика (2.23) БМЗВ внутри цилиндрического слоя (), II-асимптотика (2.24) ММЗВ (),III - гироскопическая частота ().

Вопрос классификации семейства БМЗВ - энтропийных волн намного сложен в осесимметричном случае так, как действительная часть частоты близка к нулю. Возможно, оно связано с волной БМЗВ - слоя и энтропийной модой возмущений внешней среды. Это семейство мы можем отнести к гибридным семействам БМЗВ-энтропийных волн так, как инкремент не обращается в ноль в отсутствии высвечивания (хотя в должно быть для энтропийной моды), но с другой стороны с ростом d действительная часть БМЗВ- должно стать отрицательной, все более увеличиваясь по модулю. Заметим, что в отличии от осесимметричного случая, для m=6 инкремент этого семейства гармоник оказывается крайне мал в сравнении с другими семействами.

Из выше полученных результатов, для неосесимметричных возмущений наиболее вероятно развитие гироскопических мод и гармоник семейства ММЗВ вращающегося слоя, распространяющихся в направлении его вращения. Также необходимо ответить, что инкременты отличных друг от друга гармоник любого из этих двух семейств очень слабо отличаются, поэтому неустойчивость развиваться одновременно на нескольких гармониках, другими словами, должна быть многомодовой. Последние выводы будут справедливыми в широком диапазоне параметров, если инкременты этих двух семейств слабо изменяются с изменением параметра , характеризующего интенсивность охлаждение высвечиванием, и с величиной внешнего магнитного поля (рис.2.4). При таком раскладе, наиболее крупномасштабные структуры (сравнимые с «шагом» аркады в z -направлении) должны формироваться семейством гироскопических мод.



ГЛАВА ТРЕТЬЯ. НЕУСТОИЧИВОСТЬ КОРОНАЛЬНЫХ ПЕТЕЛЬ

3.1 Модель корональной петли с продольным электрическим током

Это круглая цилиндрическая трубка радиуса a, в которой выделена центральная часть радиуса b (b < a), т. е. двойная трубка. Центральную часть назовем шнуром, а цилиндрический слой, заключенный между шнуром и внешней средой, назовем оболочкой. Плотность, давление и магнитное поле в шнуре и во внешней среде однородные, а в оболочке имеют цилиндрически-симметричные распределения:

Распределения в оболочке должны удовлетворять уравнению равновесия

Условие равновесия трубки:

Электрический ток может быть только в оболочке:

С точки зрения исследования неустойчивости интерес представляют примеры, когда имеется ненулевой продольный ток, j_0z. В общем случае линии тока есть винтовые линии, как и линии магнитного поля. Величина плотности тока меняется с радиусом, но на каждой цилиндрической поверхности r=const распределение величины тока однородное. По идеологии неустойчивость должна менять такое распределение, разрывая однородное распределение тока на множество отдельных токовых нитей. Эти нити, как ожидается, должны иметь форму, близкую к винтовым линиям. В плоском случае получаются прямые параллельные нити, лежащие на плоскости. Здесь должны получаться m винтовых нитей, расположенных на цилиндрической поверхности. Число m определяется номером цилиндрической моды.

Исследование неустойчивости целесообразно разделить на два этапа. В данной работе мы рассмотрим уравнения идеальной МГД, т. е. в отсутствие электрического сопротивления. В дальнейшем введем электрическое сопротивление.

3.2 Линейные уравнения МГД в идеально проводящей среде

Уравнения малых возмущений в идеальной покоящейся среде



где есть равновесное магнитное поле, и - равновесные давление и плотность плазмы, - возмущения поля, плотности, давления и скорости. Волны описываются в цилиндрических координатах функциями вида , где k -продольное волновое число, щ - частота, целое число m - номер цилиндрической моды.

Зависимость от радиальной переменной требует своего определения из уравнений. Радиальная компонента вектора смещения и возмущение полного давления связаны условиями [9]:

где

Далее эту систему обезразмерим, чтобы облегчить численное решение краевой задачи на ЭВМ.



Линеаризованная система МГД вне шнура (r < и r > ) имеет вид:

Данную систему свели к системе из двух обыкновенных дифференциальных уравнений вида:

где

Последнюю систему в свою очередь свели к модифицированному уравнению Бесселя:

где

Решения системы в областях 0 < r < b, b < r < a, a < r обозначим соответственно через и . Они должны подчиняться граничным условиям [10]:

3.3 Обсуждение результатов

Таким образом, выше мы поставили краевую задачу на собственные значения частоты. Еще раз выпишем ее. Получили систему из двух обыкновенных дифференциальных уравнений первого порядка для смещения и газового давления внутри шнура. Вне шнура также получили систему из двух обыкновенных дифференциальных уравнений, которая свелась к модифицированному уравнению Бесселя, что существенно облегчит решение краевой задачи. Поставленную задачу решали численно, методом стрельб на ЭВМ [12,13]. Результаты мы приводим ниже на рисунках. Результаты получены при корональных значениях безразмерных параметров задачи.

...