Лабораторная работа №5

Численное решение обыкновенных дифференциальных уравнений. Методы Рунге-Кутты.

Цель работы.

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

Краткая теория.

В практической работе №4 в методе Эйлера решения ОДУ используется член ряда Тейлора, содержащий h. В модифицированном методе Эйлера используется член ряда, содержащий h2. Ошибка этого метода составляет h3.

Математики Рунге и Кутта обобщили и развили эти методы. Согласно их классификации, метод Эйлера получил название метода Рунге-Кутты 1-го порядка. Рунге и Кутта предложили методы для решения задачи Коши, в которых более высокая точность решения дифференциальных уравнений может быть достигнута, если в методах сохранять члены ряда Тейлора, содержащие h3, h4 и т. д.

Наибольшее распространение получил метод, в котором искомая функция y(x) аппроксимируется рядом Тейлора, в котором сохранен член ряда, содержащий h4. Метод получил название метода Рунге-Кутта 4-го порядка, но в литературе он больше известен как просто метод Рунге-Кутты.

Для получения значения функции yi+1 по методу Рунге-Кутта выполняется следующая последовательность вычислительных операций:

,

где  — величина шага сетки по .

При выполнении заданий настоящей лабораторной работы мы не будем самостоятельно создавать программу, реализующую этот метод, а воспользуемся стандартной подпрограммой, использующей метод Рунге-Кутты четвертого порядка точности (подпрограмма RKGS).

Движение заряженной частицы в однородном магнитном поле.

НЕ нашли? Не то? Что вы ищете?

Рассмотрим простейший случай движения заряженной частицы в магнитном поле – плоское вращение частицы. Пусть в начальный момент времени вектор скорости электрона направлен перпендикулярно вектору индукции магнитного поля, ориентированному вдоль оси Z.

Уравнения движения электрона, записанные по компонентам скорости и координатам, в гауссовой системе единиц имеют вид:

, , , ,

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

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

Так как величина скорости в стационарном магнитном поле не меняется (здесь мы не будем рассматривать потери на излучение, тем более, что задача пока нерелятивистская) то пронормируем компоненты скорости на ее величину: , .

Учитывая связь , пронормируем координаты на .

Тогда с учётом получим безразмерный вид уравнений, описывающих движение электрона в магнитном поле:

, , , .

Решим эту систему с помощью подпрограммы RKGS.

Описание SUBROUTINE RKGS – решение системы обыкновенных дифференциальных уравнений первого порядка с заданными начальными условиями методом Рунге-Кутты 4-го порядка.

Использование в основной программе

CALL RKGS (PRMT, Y,DERY, NDIM, IHLF, FCT, OUTP, AUX) ,
где PRMT, Y,DERY, NDIM, IHLF, FCT, OUTP, AUX являются параметрами подрограммы RKGS, а параметры FCT и OUTP определяются как внешние подпрограммы.

Описание параметров.

PRMT - входной и выходной вектор длиной 5, который определяет следующие параметры:

PRMT(1)- нижний предел интегрирования системы дифференциальных уравнений.

PRMT(2)- верхний предел интегрирования системы дифференциальных уравнений.

PRMT(3)- начальный шаг интегрирования системы дифференциальных уравнений.

PRMT(4)- величина абсолютной ошибки (задаётся пользователем).

PRMT(5)- управление выводом результатов;

в начальный момент PRMT(5)=0, но если при выполнении программы PRMT(5) не равно нулю, то выполнение программы прекращается.

Y – входной вектор начальных значений функций, при выполнении программы –вектор искомых значений функций

DERY – вектор правых частей дифференциальных уравнений, но в начальный момент выполняет роль весов ошибок и в сумме должен быть равен единице.

NDIM – число уравнений в системе.

IHLF – код ошибки. Показывает, сколько раз начальный шаг интегрирования делился программой пополам (0-10), чтобы достичь заданной точности. Если IHLF больше 10 или меньше нуля, это свидетельствует о том, что заданная точность не достигается или о наличии логических ошибок в начальных условиях.

AUX – рабочий двумерный массив размером 8хNDIM.

FCT – имя подпрограммы, в которой описывается система уравнений (должна быть описана оператором EXTERNAL).Подпрограмма имеет следующие параметры: FCT(X, Y,DERY) , где Х - независимая переменная

OUTP – имя подпрограммы, в которой результаты записываются в файл выходных данных, (должна быть описана оператором external). подпрограмма имеет следующие параметры:

OUTP (X, Y,DERY, IHLF, NDIM, PRMT) , где Х - независимая переменная.

Пример использования подпрограммы RKGS приведён ниже.

! RK_Homo_B. f90

external fct, out! Описание внешних подпрограмм

real aux(8,4), u(4)

common /a/ dt, b, c, me, om, k! Общий блок памяти для основной

! программы и подпрограмм

real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, 1.0e-4, .0/)

real :: du(4) = (/0.25, 0.25, 0.25, 0.25/)

integer :: nd = 4

! Мировые константы

real :: ez = 4.8e-10, c=3.0e10, pi=3.14159, me=9.1e-28

k = 1 ! Счетчик шагов вывода результатов в файл

! Открытие файла входных данных (input. dat) и выходных данных (res. dat).

open (unit=1, file = 'input. txt',status='old')

open (unit=2, file = 'res. dat')

! Считывание входных данных из файла input. dat

read(1,*) v0, b0, al

al=al/180.*pi

dt=pi/10.

om=ez*b0/(me*c) ! Расчет частоты вращения электрона

r=v0/om ! Расчет параметра, на который нормируются координаты

v0=v0/c ! Нормировка начальной скорости электрона на скорость

! света в вакууме

b=b0/b0 ! Нормировка магнитного поля

! Безразмерное время

! Начальные величины скорости и координат

u(1)=v0*cos(al)

u(2)= v0*sin(al)

u(3)=0.0

u(4)=0.0

! Вызов подпрограммы RKGS

call rkgs(pr, u, du, nd, ih, fct, out, aux)

write(2,*) 'ih=', ih! вывод в файл результатов кода ошибки

end

! Подпрограмма, в которой описаны уравнения

subroutine fct (t, u, du)

real u(4), du(4), me

common /a/ dt, b, c, me, om, k

du(1) = - u(2)*b

du(2) = u(1)*b

du(3) = u(1)

du(4) = u(2)

return

end

! Подпрограмма, управляющая выводом результатов в файл res. dat

subroutine out (t, u, du, ih, nd, pr)

real u(4), du(4), pr(5), me

common /a/ dt, b, c, me, om, k

if (t. ge. k*dt) then

! Энергия частицы в электронвольтах

w=me*((u(1)*c)**2+(u(2)*c)**2)/2.0/1.6e-12

write (2,1) t, u(1),u(2),u(3)*(c/om),u(4)*(c/om),w!

k = k + 1

else

end if

1 format (6e12.3) ! Формат, по которому записываются результаты

! в файл res. dat

return

end

Примечание: файл входных данных input. dat создается пользователем.

Задания

1.  Протестировать программу, описанную выше, варьируя величины индукции магнитного поля, скорости электрона.

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

3.  Ввести зависимость магнитного поля от времени, задав . Провести расчеты, варьируя параметр . Объяснить полученные результаты.

4.  Ввести зависимость магнитного поля от координаты, задав . Провести расчеты, варьируя параметр . Объяснить полученные результаты.

5.  Изменить программу, обобщив её на трехмерный случай, то есть включить движение по Z. Рассчитать шаг спирали. Сравнить полученные результаты с аналитическим решением.