Поиск

Полнотекстовый поиск:
Где искать:
везде
только в названии
только в тексте
Выводить:
описание
слова в тексте
только заголовок

Рекомендуем ознакомиться

'Методические рекомендации'
Самоанализ представляет собой изучение педагогом состояния, результатов своей профессиональной деятельности, установление причинно-следственных связе...полностью>>
'Документ'
Огромная ар­мия санитаров, носильщиков, санинструкторов, фельдшеров, медсес­тер и врачей оказывала помощь раненым на переднем крае. Это персонал полк...полностью>>
'Автореферат'
Работа выполнена на кафедре источниковеденияи вспомогательных исторических дисциплин Историко-архивного институтаФедерального государственного бюджет...полностью>>
'Документ'
Ты песню, Соловей, на ветке начинайПод вечер, когда лес листом не шелохнётся:Разбужена мечтой, душа моя проснётсяИ в добрый час опять придёт весёлый ...полностью>>

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

Главная > Решение
Сохрани ссылку в одной из сетей:

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

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

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

Выполнила студентка группы 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<MAX(200,kol1) && цикл по количеству узлов (к уменьшению шага)

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

H2=h

DO WHILE tn<b && цикл по узлам отрезка времени a,b

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

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

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)<kol && цикл по количеству уравнений

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

k1&j=f(j,1)

ENDDO

j='0'

DO WHILE val(j)<kol

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

k2&j=f(j,2)

ENDDO

j='0'

DO WHILE val(j)<kol

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

k3&j=f(j,3)

ENDDO

j='0'

DO WHILE val(j)<kol

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.



Скачать документ

Похожие документы:

  1. М. В. Ломоносова Факультет вычислительной математики и кибернетики В. Г. Баула Введение в архитектуру ЭВМ и системы программирования Москва 2003 Предисловие Данная книга

    Книга
    Данная книга представляет собой учебное пособие по архитектуре ЭВМ и системам программирования. Книга написана по читаемому автором лекционому курсу Архитектура ЭВМ и язык Ассемблера для студентов первого курса факультета Вычислительной
  2. М. В. Ломоносова Факультет вычислительной математики и кибернетики Руденко Т. В. Сборник задач

    Сборник задач
    Представлены задачи и упражнения по языку Си и программированию на нем. Рассматриваемая версия Си соответствует международному и ANSI-стандарту этого языка.
  3. М. В. Ломоносова Факультет вычислительной математики и кибернетики Н. В. Вдовикина, А. В. Казунин, И. В. Машечкин, А. Н. Терехин Системное программное обеспечение: взаимодействие процессов учебно-методическое пособие

    Учебно-методическое пособие
    В пособии рассматриваются основные аспекты управления процессами в операционной системе и организации межпроцессного взаимодействия на примере операционной системы UNIX.
  4. И кибернетики факультет вычислительной математики и кибернетики

    Доклад
    Доклад доцента Панкратова А.Н., доцента Поповой Н.Н., профессора Дедуса Ф.Ф., ст.науч.сотр. ИМПБ РАН Тетуева Р.К., асп.Пяткова М.И., асп. Колесина М.С.
  5. Московский государственный университет имени М. В. Ломоносова Факультет вычислительной математики и кибернетики

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

Другие похожие документы..