Geum.ru

Факультет вычислительной математики и кибернетики

Вид материалаРешение

Содержание


2.Метод решения
3.Интерфейс программы
4.Выполнение программы
5.Тест программы
6.График орбиты Аренсторфа
7.Основной текст программы (P1.prg)
Подобный материал:


КАЗАНСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

ФАКУЛЬТЕТ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И КИБЕРНЕТИКИ


Семестровая работа по численным методам: «Решение системы дифференциальных уравнений методом Рунге-Кутта. Ограниченная задача трех тел»


Выполнила студентка группы 916

Габидинова Айгуль Ринатовна.

Преподаватель: Задворнов О.А.


Казань-2004

1.Постановка задачи


Рассмотрим два тела с массами m (Луна) и M=1-m (Земля), участвующие в совместном круговом движении в плоскости xOy и расположенные в точках с координатами (1,0)
и (0,0) соответственно. Пусть далее вблизи этих тел в той же плоскости движется третье тело пренебрежимо малой массы и (x(t),y(t)) его координаты в момент времени t. Траектория движения этого тела описывается уравнениями:


x″=x+2y′-M(x+m)/R1-m(x-M)/R2

y″=y-2x′-My/R1-my/R2 (1)

R1=((x+m)2+y2)3/2, R2=((x-M)2+y2)3/2,


Уравнения (1) дополняются начальными условиями:

x(0)=0.994, y(0)=0, x′(0)=0, y′(0)=-2.031732629557337. (2)


При начальных условиях (2) и m=0.012277471 орбита будет периодической с периодом обращения равным T=11,124340337 (такие орбиты называют “орбитами Аренсторфа”).

2.Метод решения


Для начала сведем заданные уравнения к виду системы из 4 уравнений первого порядка, обозначив p=x′, q=y′:


x′=p

y′=q (5)

p′=x+2q-(1-m)(x+m)/((x+m)2+y)3/2-m(x-1+m)/((x-1+m)2+y)3/2

q′=y-2p-(1-m)y/((x+m)2+y2)3/2-my/((x-1+m)2+y2)3/2

с начальными условиями

x(0)=0.994, y(0)=0,p(0)=0, q(0)=-2.031732629557337


Система из n дифференциальных уравнений на произвольном отрезке [0,T] решается методом Рунге-Кутта четвертого порядка точности:

k1=f(tn,yn), k2=f(tn+h/2,yn+hk1/2),

k3=f(tn+h/2,yn+hk2/2), k4=f(tn+h,yn+hk3), (3)

yn+1=yn+h(k1+2k2+2k3+k4)/6.

Программа написана в СУБД VISUAL FOXPRO8, проект и EXE файл называются KOSHI.

С экрана запрашиваются начало A и конец B отрезка интегрирования и количество узлов KOL1. Шаг интегрирования вычисляется в программе H =(B-A)/KOL1


Процедура интегрирования в программе P1 по формулам (3) реализуется следующим образом: сначала вычисляются коэффициенты k1 поочередно для всех решений системы уравнений в цикле по количеству уравнений в системе, в данном случае 4 (для тестовой задачи 2), затем аналогично k2, k3 и k4. Формулы уравнений реализуются процедурой-функцией F для x, y, p, q по формулам (5) (для тестовой задачи для y1, y2 по формулам (4)). Входными параметрами этой функции F(j,d) являются j – порядковый номер решения системы уравнений и d- порядковый номер коэффициента.


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

y1′=-y2+ y1 (y12+y22-1)

y2′=y1+ y2 (y12+y22-1) (4)

на отрезке [0,5] с точным решением

y1=cos(x)/(1+e2x)1/2 , y2=sin(x)/(1+e2x)1/2

В результате получаем max погрешность порядка 10-5.

3.Интерфейс программы


Главная форма программы состоит из экранной формы с тремя вкладками: “Тестовая задача Коши”, “Ограниченная задача трех тел” и “Орбита”, в двух присутствуют запросы отрезка интегрирования, количества узлов и кнопка “ОК”, по нажатию которой осуществляется запуск соответствующего варианта программы.

4.Выполнение программы


В варианте тестовой задачи программа запускается в цикле от 25 до введенного количества узлов, рассчитываются соответствующие шаги H, в каждом случае вычисляются точное и приближенное решения, вычисляется max погрешность по всем временным слоям max(y1точное -y1приближенное, y2точное -y2приближенное). Данные по шагам и соответствующим погрешностям заносятся в файлы POGR.DBF и TEST.XLS, а вычисленные решения при заданном с экрана количестве узлов заносятся в файлы REZULT.DBF и N.XLS.

В варианте ограниченной задачи трех тел программа запускается один раз для заданного с экрана количества узлов, рассчитанные координаты заносятся в файлы REZULT.DBF и ARENSTORF.XLS. В окрестностях начала и конца периода выбираем существенно меньший шаг интегрирования H, чем в другие моменты времени.

Файлы DBF выводятся на экран в табличном виде по ходу работы программы. Файлы XLS позволяют получить в EXCEL графики.

5.Тест программы


Для проверки правильности работы процедуры интегрирования задачи Коши методом Рунге-Кутта была решена тестовая задача.

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







6.График орбиты Аренсторфа


Вид траектории движения тела в системе координат (x,y) относительно Земли и Луны на отрезке интегрирования [3,26] и количестве узлов 4600:


7.Основной текст программы (P1.prg)


SET PROCEDURE TO p1

CREATE dbf rezult (t n(18,16),y1 n(20,16),y2 n(20,16),y3 n(20,16),y4 n(20,16),z1 n(20,16),z2 n(20,16),pogresh n(20,16))

CREATE dbf pogr (shag n(18,16),pogresh n(20,16),pogr_ n(20,16))

SELECT 2

USE pogr

SELECT 1

USE rezult

zad=1 && zad= 1 - тестовая задача, 2 - ограниченная задача трех тел

a=0 && a,b - отрезок времени

b=5

c=11.124340337

period=c

kol1=100 && при kol1 количестве узлов выдать таблицу результатов

kol2=0

DO FORM form1.scx

kol=IIF(zad=1,2,4) && система из kol уравнений

IF zad=2

m=0.012277471

b=c

kol1=kol2

ENDIF

kol1=MAX(200,kol1)

N=IIF(zad=1,25,kol1-1) && количество узлов

DO WHILE N
SELECT 1

ZAP

N=N+1

tn=a

h=(b-a)/N && шаг для данного количества узлов N

APPEND BLANK && заполним начальные данные

IF zad=1

REPLACE t with tn,y1 with cos(a)/sqrt(1+exp(2*a)),y2 with sin(a)/sqrt(1+exp(2*a))

ELSE

REPLACE t WITH tn,y1 WITH 0.994,y2 WITH 0,y3 WITH 0,y4 WITH -2.031732629557337

ENDIF

h1=h

DO WHILE tn
h=IIF(zad=2.and. (mod(t-a,period)<0.005*period.or.mod(t-a,period)>period*0.995),0.00001,h1) && выбираем существенно меньший шаг в окрестностях точек, кратных периоду

*h=IIF(zad=2.and. (mod(t-a,period)<2*h1.or.mod(t-a,period)>=period-h1),0.00005,h1) && выбираем существенно меньший шаг в окрестностях точек, кратных периоду

tn=tn+h && очередной узел

u1=y1 && запомним предыдущие y1,y2,y3,y4 для вычисления следующих

u2=y2

u3=y3

u4=y4

APPEND BLANK

REPLACE t WITH tn

* вычисления методом Рунге-Кутта 4 порядка

j='0'

DO WHILE val(j)
j=str(val(j)+1,1)

k1&j=f(j,1)

ENDDO

j='0'

DO WHILE val(j)
j=str(val(j)+1,1)

k2&j=f(j,2)

ENDDO

j='0'

DO WHILE val(j)
j=str(val(j)+1,1)

k3&j=f(j,3)

ENDDO

j='0'

DO WHILE val(j)
j=str(val(j)+1,1)

k4&j=f(j,4)

REPLACE y&j WITH u&j+h*(k1&j+2*k2&j+2*k3&j+k4&j)/6

ENDDO

ENDDO

IF zad=1

GO TOP

max_=0

* вычисление точного значения (z)

DO WHILE .not. eof()

REPLACE z1 with cos(t)/sqrt(1+exp(2*t))

REPLACE z2 with sin(t)/sqrt(1+exp(2*t))

REPLACE pogresh WITH max(abs(y1-z1),abs(y2-z2)) && вычисление погрешности sqrt((y1-z1)**2+(y2-z2)**2)

max_=max(max_,pogresh) && фиксируем max погрешность для данного шага

SKIP

ENDDO

IF n=kol1 && вывод на экран таблицы результатов при заданном количестве узлов kol1

COPY TO n TYPE XL5 FIELDS pogresh

BROWSE TITLE 'Результаты при количестве '+STR(kol1,5)+' (n.xls)'fiel t :h='Узлы',pogresh :h='Погрешность',y1 :h='y1 прибл',;

z1 :h='y1 точн',y2 :h='y2 прибл',z2 :h='y2 точн' STYLE 'B' NODELETE NOMENU

ENDIF

SELECT 2

APPEND BLANK && заполнение таблицы шагов и погрешностей

REPLACE shag with h,pogresh with max_,pogr_ with pogresh/h**4

ELSE && zad=2

COPY TO arenstorf TYPE XL5 FIELDS t,y1,y2,y3,y4

BROWSE TITLE 'Координаты орбиты X,Y при количестве узлов'+STR(kol1,4)+'(k.xls)' fiel t :h='Узлы',;

y1 :h='X',y2 :h='Y' STYLE 'B' NODELETE NOMENU

ENDIF zad=1

ENDDO

IF zad=1

SELECT 2

COPY TO test TYPE XL5 && копирование таблицы в M.XLS в формате EXCEL для построения графиков

* вывод на экран таблицы зависимости шагов и погрешностей

BROWSE TITLE 'Зависимость MAX погрешности от шага (m.xls)' fiel shag :h='Шаг',pogresh :h='max погрешность',;

pogr_ :h='max погрешность/шаг в 4 степ' STYLE 'B' NODELETE NOMENU

ENDIF

RETURN

*процедуры формул производных, по которым вычисляются решения методом Рунге-Кутта

PROCEDURE f

PARAMETERS j_,d_

DO CASE

CASE j_='1' .and.zad=1

RETURN -(u2+d(2))+(u1+d(1))*((u1+d(1))**2+(u2+d(2))**2-1) && для тестового

CASE j_='2'.and.zad=1

RETURN (u1+d(1))+(u2+d(2))*((u2+d(2))**2+(u1+d(1))**2-1) && для тестового

* для задачи трех тел

CASE j_='1'.and.zad=2

RETURN u3+d(3) && для огран.задачи трех тел

CASE j='2'.and.zad=2

RETURN u4+d(4) && для огран.задачи трех тел

CASE j_='3'

RETURN (u1+d(1))+2*(u4+d(4))-(1-m)*(u1+d(1)+m)/SQRT(((u1+d(1)+m)**2+(u2+d(2))**2)**3)-m*(u1+d(1)-(1-m))/SQRT(((u1+d(1)-1+m)**2+(u2+d(2))**2)**3)

CASE j_='4'

RETURN (u2+d(2))-2*(u3+d(3))-(1-m)*(u2+d(2) )/SQRT(((u1+d(1)+m)**2+(u2+d(2))**2)**3)-m*(u2+d(2) )/SQRT(((u1+d(1)-1+m)**2+(u2+d(2))**2)**3)

ENDCASE

PROCEDURE d

PARAMETERS l

l=STR(l,1)

RETURN IIF(d_<3,IIF(d_=1,0,h*k1&l/2),IIF(d_=3,h*k2&l/2,h*k3&l))


Литература:

  1. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. – М.: Мир, 1990.
  2. Самарский А.А., Гулин А.В. Численные методы. – М.:Наука, 1989.
  3. Базвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. – М.: Наука, 1987.