СКОРОСТЬ ГРАВИТАЦИИ

             Часть 3 – АНОМАЛЬНЫЕ СМЕЩЕНИЯ ПАРАМЕТРОВ ОРБИТ ПЛАНЕТ

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

                                                      ….первая редакция 29.12.2012                                                                                  

 

С. Ю. Юдин................................  …http://modsys.narod.ru.. ................. ..modsys@narod.ru

 

 

В середине 19-го века Леверье при создание своей теории планет обнаружил, что наблюдаемые значения положений планет отличаются от рассчитанных им по теории Ньютона. Но, если в этой теории планет учесть систематические поправки по смещению перигелиев Меркурия и Марса, соответственно, на 38,3 и 25,15 угловых секунд, то это расхождение устраняется. В 1895 году Ньюком при создание своей теории планет, уточняя таблицы Леверье, тоже обнаружил эти расхождения, а также обнаружил расхождения и по положению перигелиев других планет и других параметров орбит планет. Но достоверными он признал только расхождения, т.е. аномальные отклонения, в вековых смещениях перигелиев Меркурия и Марса, узлов Венеры и эксцентриситета Меркурия. В дальнейшем эти величины многократно в различных теориях уточнялись, но, например, французы до сих пор пользуются теориями созданными их соотечественником Леверье. При этом, начиная с Ньюкома, астрономы стараются не просто вносить эмпирические поправки к расчетным значениям параметров орбит планет, а сделать так, чтобы эти поправки вытекали из примененных ими при создании математической модели движения планет физических законов.

 

 

Так Ньюком использовал при этом модифицированный (с поправкой Холла) закон Ньютона, который давал дополнительное смещение перигелиев планет. И хотя Холл, предложивший определять поправку для закона Ньютона (1) именно по самому заметному смещению, т.е. смещению перигелия Меркурия, и нашел значение этой поправки n=0,0000001552, Ньюком принял другое значение этой поправки n=0,0000001612. В связи с этим по его теории аномальные смещения теперь отличаются от им самим же полученных аномальных смещений перигелиев планет (см. табл. 1). В этой же таблице Вы видите и результаты по некоторым другим теориям и не только физическим (теория Зеелингера чисто астрономическая), которые тоже позволяют сместить перигелии планет на угол больше, чем получается по теории Ньютона. И хотя, как показал Ньюком, аномальные смещения имеются и у других параметров орбит, но ученые того времени сосредоточились в основном на объяснение аномального смещения перигелия Меркурия, а тут, как мы знаем, преуспела ОТО Эйнштейна.

 

F = gamma * m0 * m / R^(2+n)      (1)

 

Хотя, как мы видим из таблицы, все представленные в ней теории отлично объясняют аномальное смещение перигелия Меркурия и гораздо хуже смещение перигелиев других планет. Но все они (кроме теории Зееленгера) никак не объясняет аномальные смещения линии узлов, углов наклона орбит и их эксцентриситетов, которые по данным Ньюкома (табл. 10) тоже имеются, по этому, если какая то теория объяснит не только аномальное смещение перигелия Меркурия, но и аномальные смещения других параметров орбит планет, то эта теория будет более истинна, чем принятая сейчас за эталон истины ОТО. Правда в работе [14] говорится, что и ОТО дает некоторое дополнительное смещение перигелиев и узлов планет от вращения Солнца, но это смещение столь мало, что никакого практического влияния на расчетное смещение перигелиев и узлов планет оказать не может. Например, по данным  [14] смещение перигелия Меркурия от вращения Солнца будет -0,01 угл. сек., а смещение узлов будет еще в 2 раза меньше, что практически никак не повлияет на аномальные остатки в таблице 10. Таким образом, можно считать, что и ОТО никак не объясняет аномального смещения узлов планет, углов наклона и эксцентриситетов их орбит. Следовательно, надо дальше уточнять наблюдательные данные и создавать теории, которые объяснят все наблюдаемые смещения параметров орбит.

 

 

Таблица 1. Аномальное смещение перигелия внутренних планет (в угл.сек. за 100 лет).

 

 

Меркурий

Венера

Земля

Марс

Наблюдения  (Ньюком [6])

575,06

42,52

1162,92

1602,69

теория Ньютона (Ньюком [6])

533,82

49,85

1156,95

1594,65

Аномальный остаток (Ньюком [6])

41,24

-7,33

5,97

8,04

теория Холла (Ньюком [7])

43,37

16,98

10,45

5,55

теория Эйнштейна (Субботин [7])

43,03

8,62

3,83

1,35

теория Гербера (Хайдаров [8])

43,03

8,62

3,83

1,35

теория Ритца (Роузвер [9])

41,0

8,0

3,4

-

теория Маха (Зайцев [10])

43,0

23,0

17,0

11,0

теория Зеелингера (Роузвер [9])

41,3

7,3

4,2

6,3

 

 

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

 

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

 

Таким образом, в 20-м веке появилось множество возможностей для уточнения аномальных смещений параметров орбит планет, которые можно подкорректировать не только уточнением наблюдаемых смещений параметров орбит, но и расчетных, т.е. посчитанных по теории Ньютона путем изменения некоторых астрономических постоянных. И ниже я привожу в табл. 2 некоторые системы масс планет, которые использовались при создании тех или иных теорий планет для расчета смещений параметров орбит по теории Ньютона. А в табл.3, я показываю как может измениться расчетное значение смещения перигелия Меркурия при использование различных систем масс планет. Я рассчитывал смещение на компьютере и мог сразу посчитать суммарное значение, а не вклад каждой планеты, как это был вунужден делать Леверье решая аналитически задачу двух тел, но для наглядности я сделал так как было у него (смещения определялись по методике с критерием максимина за период с 1601 по 2000 год). При этом, естественно, Леверье использовал свою систему масс, а я использовал систему масс JPL и, как я догадываюсь, Клеменс использовал систему масс Ньюкома.  При этом, даже современную систему масс JPL следует признать не идеальной, т.к., как мы видим в табл.3, наибольшее влияние на смещение перигелия Меркурия оказывают Венера, Земля и Юпитер, но до сих пор их массы определены не точно. Согласно данным, приведенным в [11], у Венеры масса находится в интервале от 408 250 (Эш, Шапиро, Смит – 1967) до 408 945 (Данком – 1965), у Земли от 328 890 (Рабе – 1967) до 328 939 (Рабе, Френсис – 1967) и у Юпитера от 1047,314 (Черных – 1970) до 1047,445 (Клепчинский – 1964). И это не смотря на то, что и у Земли и у Юпитера есть спутники. Что уж говорить о Венере, у которой спутника нет.

 

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

 

Леверье

Ньюком

МАС 1964

JPL

Меркурий

3 000 000

6 000 000

6 000 000

5 983 000

Венера

401 847

408 000

408 000

408 522

Земля+Луна

354 936

329 390

329 390

328 900,1

Марс

2 680 337

3 093 500

3 093 500

3 098 700

Юпитер

1 050

1 047,355

1 047,355

1 047,391

Сатурн

3 512

3 501,6

3 501,6

3 499,2

Уран

24 000

22 869

22 869

22 930

Нептун

14 400

19 314

19 314

19 600

Плутон

-

-

360 000

1 812 000

 

 

Таблица 3. Расчетные значения смещения перигелия Меркурия от действия отдельных планет с разными системами масс (в угл.сек. за 100 лет).

 

 

Леверье 1856 [3, 7]

Клеменс 1947 [1]

Юдин 2008

Меркурий*

-

+0,025

+0,0032

Венера

+280,64

+277,856

+275,95

Земля

+83,61

+90,038

+90,27

Марс

+2,85

+2,536

+2,48

Юпитер

+152,59

+153,584

+153,13

Сатурн

+7,24

+7,302

+7,24

Уран

+0,14

+0,141

+0,14

Нептун

+0,06

+0,042

+0,05

сжатие Солнца

-

+0,01

-

итого

+527,13

+531,534

+529,29


* - у меня это ошибка от численного решения дифференциальных уравнений и погрешности при статистической обработке данных, т.е. смещение перигелия Меркурия когда он один вращается вокруг Солнца, а вот что это у Клеменса я затрудняюсь сказать, т.к. в 1947 году, когда по данным Брумберга были опубликованы эти данные, хоть и появились допотопные ЭВМ, но на них навряд ли можно было произвести такие сложные расчеты для Меркурия. Хотя возможно, что и это тоже ошибка численного решения, т.к. уже в 1952 году Брауэр и Клеменс в Морской обсерватории США  решали на ЭВМ численными методами систему уравнений для внешних планет (правда, они  делали это с очень большим шагом - 40 дней, чтобы существующие уже в 1952 году ЭВМ хоть как то справились с этой задачей).

 

Такая же неопределенность существует на сегодняшний день и с наблюдаемыми смещениями параметров орбит. Но если полученные тем или иным автором расчетные значения вековых смещений параметров орбит при наличие компьютера легко проверить, то с наблюдаемыми данными все гораздо сложнее. И даже, когда Дулитл в 1918 году решил перепроверить расчетные значения смещения перигелия Меркурия у Леверье и Ньюкома он столкнулся с большой неожиданностью. Дело в том, что Леверье и Ньюком определяли расчетные значения по менее точному методу Лапласа, а Дулитл использовал более точный метод Гаусса по которому у него для Меркурия получалось 529,67 угл.сек. Тогда он повторил расчеты Леверье и Ньюкома тоже по методу Гаусса и с принятой им (Дулитлом) системой масс планет и получил, соответственно, 529,84 и 533,03 угл.сек. [3] Расхождение с данными Леверье составило 0,17, что вполне допустимо, а вот с данными Ньюкома расхождение получилось 3,36. В чем тут причина Дулитл так и не понял, но в 1926 году Шази указал [3], что это получилось из-за того, что Ньюком пользовался не правильной формулой для определения смещения перигелия, и по этому у него в расчетах имеется ошибка как раз на 3,36 угл.сек. и, следовательно, расчетное значение у него должно было быть 530,46, а не 533,82, как у нас указано в табл. 1.     

 

 

Учитывая такой разброс даже в расчетных данных, я ниже приведу некоторые данные и по наблюдаемым и по расчетным значениям вековых смещений перигелиев 4-х внутренних планет, как с расчетными данными самих авторов, так и со своими расчетными данными, которые дам в знаменателе. Свои расчетные данные я определял сразу от действия всех планет, а не по отдельности, как я это делал в табл.3. При этом, например, данные, полученные мною суммарно для системы масс JPL (529,67), немного отличаются от раздельного расчета в табл.3 (529,29), но это обусловлено погрешностями статистической обработки расчетных данных по методике которую я изложу ниже. В расчете для данных JPL значение получено с системой масс JPL, а для остальных данных с системой масс Ньюкома. При этом данные Ньюкома я приведу в двух вариантах – так как их привел сам автор и полученные мною по современной версии его теории, которая опубликована в [11].

 

Таблица 4-1. Наблюдательные и расчетные данные по смещению перигелия Меркурия по расчетным данным самих авторов в числителе и по моим расчетным данным в знаменателе.

 

наблюдения

расчет

остаток

Ньюком [6]

575,06

533,82 / 529,81

41,24 / 45,25

Ньюком [11]

570,73

- / 529,81

- / 40,92

Данком [1]

575,15

532,05 / 529,81

43,10 / 45,34

JPL0 [2]

577,73

- / 529,67

- / 48,06

JPL2 (DE405) [4]

572,20

- / 529,67

- / 42,53

 


Таблица 4-2. Наблюдательные и расчетные данные по смещению перигелия Венеры по расчетным данным самих авторов в числителе и по моим расчетным данным в знаменателе.

 

наблюдения

расчет

остаток

Ньюком [6]

42,52

49,85 / 45,33

-7,33 / -2,81

Ньюком [11]

54,41

- / 45,33

- / 9,08

Данком [1]

34,29

26,22 / 45,33

8,07 / -11,04

JPL0 [2]

9,66

- / 44,79

- / -35,13

JPL2 (DE405) [4]

43,64

- / 44,79

- / -1,15

 


Таблица 4-3. Наблюдательные и расчетные данные по смещению перигелия Земли по расчетным данным самих авторов в числителе и по моим расчетным данным в знаменателе.

 

наблюдения

расчет

остаток

Ньюком [6]

1162,92

1156,95 / 1144,02

5,97 / 18,90

Ньюком [11]

1160,37

- / 1144,02

- / 16,35

Данком [1]

1153,77

1148,76 / 1144,02

5,01 / 9,75

JPL0 [2]

1163,77

- / 1142,66

- / 21,11

JPL2 (DE405) [4]

1158,55

- / 1142,66

- / 15,89

 

 

Таблица 4-4. Наблюдательные и расчетные данные по смещению перигелия Марса по расчетным данным самих авторов в числителе и по моим расчетным данным в знаменателе.

 

наблюдения

расчет

остаток

Ньюком [6]

1602,69

1594,65 / 1602,17

8,04 / 0,52

Ньюком [11]

1601,09

- / 1602,17

- / -1,08

JPL0 [2]

1599,88

- / 1598,00

- / 1,88

JPL2 (DE405) [4]

1600,02

- / 1598,00

- / 2,02

 

 

Как мы видим по данным таблиц 4-1...4-4 для аномального смещения перигелиев планет у нас имеется широкий выбор данных. Тоже самое мы имеем и по аномальным смещениям других параметров орбит планет. Например, как я писал выше, Ньюком в 1895 году обнаружил аномалию в смещение узлов Венеры, а Данком в 1958 году ее не обнаружил. А если к этому добавить влияние на результат не только разных систем масс, но и других астрономических постоянных, например, величины астрономической единицы или влияние сжатия Солнца, то результаты получатся еще более пестрые. Ведь, например, согласно данным 1967 года Дикке и Голденберга [1, 9] сжатие Солнца увеличивает вековое смещение перигелия Меркурия на 3,6 угл.сек., а во всех приведенных выше данных (кроме Клеменса в табл.3) эта поправка не учитывалась. Вот по этому я в предыдущей статье [4] и создавал свою кинематическую теорию планет, чтобы найти чисто наблюдательные, т.е. не завязанные ни на какую теорию, вековые смещения параметров орбит по которым я буду искать скорость распространения гравитации и скорость движения Солнечной системы. И ниже я привожу табл. 5 с найденными мною значениями этих величин (таблица приводилась в [4]). Смещения углов здесь в угл.сек., а смещения эксцентриситетов я привожу не в угловых секундах, как это делал Ньюком, а в безразмерных величинах увеличенных в 1 000 000 раз.

 

Таблица 5. Вековые смещения параметров орбит, полученные на программах Solsys6 и Solsys7 с использованием данных полученных по имитаторам (астрономическим теориям) Солнечной системы New0, APC2 и JPL2 и данным математических моделей Солнечной системы (в интервале с 1601 по 2000 год), а также полученных при обработке непосредственно данных оптических наблюдений для создания кинематической астрономической теории Ser0 в интервале примерно с 1800 по 1980 год.

 

 

параметры

имитаторы

модели

New0

APC2

JPL2

Ser0

/R^2

/R^n*

dAlfaP1

dAlfaU1

dBetta1

dEks1

+570,73

-452,18

-21,43

+20,55

+574,02

-450,34

-21,44

+20,52

+572,20

-449,95

-21,44

+20,50

+578,0

-433,2

-19,84

+20,10

+529,59

-450,36

-21,455

+20,421

+573,41

-450,40

-21,457

+20,417

dAlfaP2

dAlfaU2

dBetta2

dEks2

+54,41

-999,86

-2,51

-48,01

+59,12

-1011,65

-2,42

-48,20

+43,64

-998,36

-2,52

-48,43

30,3           -996,5         -2,83         -45,33

+45,555

-999,08

-2,521

-49,417

+61,193

-999,06

-2,521

-49,409

dAlfaP3

dAlfaU3

dBetta3

dEks3

+1160,37

-869,85

-47,19

-41,33

+1185,68

-872,74

-47,14

-41,62

+1158,55

-861,63

-47,17

-41,93

+1141,6

-868,0

-47,45

-41,57

+1141,97

-845,52

-47,235

-43,291

+1152,56

-845,36

-47,234

-43,291

dAlfaP4

dAlfaU4

dBetta4

dEks4

+1601,09

-1054,74

-29,11

+92,27

+1599,64

-1058,58

-29,16

+92,07

+1600,02

-1050,49

-28,93

+91,63

1582,7       -1028,0     -29,39        +95,53

+1598,09

-1052,79

-29,025

+96,255

+1603,70

-1052,77

-29,025

+96,252

 

* - показатель степени n во второй модели в формуле Ньютона F=G*m*M/R^n брался таким, каким он был у Ньюкома, т.е. n=2,0000001612 (у Холла было n=2,0000001574).

 

 

Но перед тем, как приступить к поиску скорости гравитации, надо выяснить какие параметры орбит планет более всего отклоняются от параметров, которые получаются при моделирование Солнечной системы по закону притяжения Ньютона, где скорость гравитации принимается равной бесконечности. Для этого надо не только найти аномальные остатки, но и выяснить какова достоверность этих отклонений, т.е. с какой точностью мы их определили, а эта точность будет зависеть от точности определения наблюдаемых и расчетных смещений параметров орбит. И вот тем как уменьшить доверительные интервалы хотя бы для расчетных данных мы сейчас и займемся. Для определения расчетных значений параметров орбит можно бы было использовать ту же методику, что я использовал при обработке данных наблюдений в предыдущей статье [4], но это очень трудоемкий процесс. А учитывая то, что расчетных данных у нас при поиске скорости гравитации будет не просто много, а очень много, данная методика меня не устраивает. По этому я для обработки расчетных данных, т.е. просто координат X, Y, Z, которые мы получаем при моделирование движения планет в Солнечной системе с использованием различных физических законов, разработал такую двухступенчатую методику. На форме 2, где моделируется движение планет мы отмечаем не чекбокс //Записать координаты X, Y, Z//, чтобы потом обрабатывать их на фоме 10, как данные оптических наблюдений, и только потом переходить к форме 6, где производится статистическая обработка данных по полученным параметрам орбит, а чекбокс //Записать параметры орбит в файл//. При этом параметры орбит у нас определяются непосредственно во время моделирования при каждом обороте планеты и потом уже этот массив данных мы будем обрабатывать на форме 6.

 

 

Рис. 1. Моделирование движения планет Солнечной системы и определение мгновенных параметров орбиты Меркурия. Скриншот программы Solsys7.

 

Таким образом, первая часть этой методики состоит в том, чтобы определить во время движения планеты параметры орбиты, по которой она должна так двигаться. Угол восходящего узла AlfaU определяется просто. Как только координата Z стала не отрицательной, а положительной мы присваиваем AlfaU текущее значение эклиптической долготы планеты, которое легко определяется по текущим координатам X и Y. Несложно во время движения определить и среднюю угловую скорость Wsr и угол наклона орбиты Betta. Wsr мы находим из периода обращения планеты, а его находим по двум моментам времени, когда эклиптическая долгота равна нулю сейчас и когда она была равна нулю при предыдущем обороте планеты. А угол наклона орбиты принимаем тот который был максимальным из всех Betta определяемых на каждом шаге решения дифференциальных уравнений при одном обороте планеты. Немного сложнее с параметрами Eks (эксцентриситет) и R (большая полуось). Здесь нам надо найти максимальный и минимальный радиусы орбиты, которые были при одном обороте планеты, а затем уже по этим значениям рассчитать Eks и R. Сложнее всего определить мгновенное значение перигелия AlfaP. Здесь нам надо не просто определить минимальное значение R, а зафиксировать момент, когда R после того как он уменьшался начал увеличиваться и в этот момент мы присваиваем AlfaP текущее значение угла равного сумме двух углов. Первый угол это, определенный ранее в плоскости эклиптики, угол восходящего узла AlfaU, т.е. угол между осью координат Х и линией пересечения плоскости орбиты с плоскостью эклиптики, а второй угол это угол в плоскости движения планеты от этой линией пересечения до положения планеты на своей орбите, которая наклонена к плоскости эклиптики под найденным ранее углом Betta.

 

Астрономы могут мне возразить, что это мгновенные параметры орбиты, т.е. параметры оскулирующей орбиты, которые постоянно изменяются и по этому они будут отличаться от невозмущенной орбиты планеты. Согласен, и на рис. 1, где отмечен чекбокс dAlfaP, мы видим, что на график выводится смещение перигелия Меркурия (ломаная линия из белых точек от начала оси абсцис), которая не является идеальной прямой. Но ведь в общем случае мы видим, что это смещение увеличивается линейно, а мгновенные значения параметра просто отклоняются как случайные величины от этой глобальной прямой. И вот эту глобальную зависимость нам и надо найти. Для этого мы можем обработать методом наименьших квадратов эти данные на форме 6, как мы это делали при обработке данных наблюдений полученных на форме 10, и получить глобальную зависимость и предельные отклонения от этой зависимости, но мы будем поступать несколько иначе. Давайте рассмотрим тестовый пример на форме 6 (рис. 2), где у нас 25 значений параметра рассчитываются по формуле (2), если отмечен чекбокс у этой формулы.

 

Alfa = k0 + k1 * dT + k2 * dT^2 + k3 * Rnd (dT – 0,5)                  (2)

 

 

Рис. 2. Статистическая обработка данных тестового примера на форме 6. Скриншот программы Solsys 7.

 

На рисунке Вы видите 25 красных точек рассчитанных по формуле (2), когда k0=0, k1=1, k2=0, k3=1 и все эти точки обведены большим черным кружком, который означает, что данные обрабатывались в исходном виде. Теперь разобьем всю выборку на 5-ть групп, вычислим средние значения параметра по 5-и точкам в группе и потом обработаем эти средние 5-ть точек (см. рис. 3). Теперь большими черными кружками отмечены средние значения параметра в группах. Коэффициенты уравнения регрессии получились примерно одинаковыми в обоих случаях, а вот дисперсия во втором случае получилась более чем в два раза меньше. Причем коэффициенты аппроксимации в двух случаях просто не могут получиться одинаковыми по той простой причине, что коэффициент k3 в уравнение (2) умножается на случайное число, которое генерирует функция Rnd. При этом, хотя дисперсия и получается немного разной, но она устойчиво, более чем в два раза меньше во втором варианте аппроксимации данных.

 

 

Рис. 3. Статистическая обработка данных тестового примера на форме 6. Скриншот программы Solsys 7.

 

Строгого математического обоснования правомочности такой обработки данных по средним значениям в группах я не нашел, но проведенный нами эксперимент говорит о том, что результат всегда получается точнее, т.е. с меньшим доверительным интервалом при одинаковой доверительной вероятности (надежности). По этому, я буду использовать этот эффект при обработке данных по параметрам орбит полученным при моделирование движения планет. А теперь давайте обработаем с учетом этого эффекта данные по расчетным параметрам орбит, которые будем определять не по результатам работы математической модели, а по данным теории JPL2, т.е. по данным, которые будем рассчитывать по эфемеридам DE405, которые в свою очередь являются очень точной аппроксимацией (отрезками по 8...32 дня) полиномами Чебышева результататов работы математической модели JPL, т.е. по таким же данным, которые бы мы получили отметив чекбокс //Записать координаты X, Y, Z// при работе нашей модели. Но на форме 2 теория JPL2 отнесена к группе имитаторов, т.е. симуляторов или если хотите симулянтов выдающих себя за модели. Такими же имитаторами являются и теории New0, JPL0 и Ser0 описанные мною в предыдущей статье. Более подробно о классификации моделей и имитаторов читайте в моих работах [5, 11].

 

Обработать данные с формы 2 мы можем или разбив всю первичную выборку данных на группы и обработав вторичную выборку как средние значения в группах или как отдельные наблюдения, и при этом второй расчет будет одним из вариантов разбивки данных на группы, когда в группе одно наблюдение. Возникает только вопрос на какое количество групп разбить первичную выборку данных, чтобы получить наиболее достоверный результат. На рис. 4 Вы видите матожидания векового смещения параметра AlfaP (синяя кривая) и верхнее (красная кривая) и нижнее (зеленая кривая) предельные отклонения этого параметра от матожидания, которые определялись как 2 стандарта, при разном количестве точек в группе от 1 (начало) до 1106 (конец) при масштабе по оси абсцисс 50 точек в сантиметре (в первичной выборке 3318 значений). Как видим при малом количестве точек в группе у нас получается большая дисперсия, которая при увеличение точек в группе не просто уменьшается, а еще и периодически убывает и возрастает.  Это связано с тем, что при разном количестве точек в группе, т.е. при разном периоде времени, за который мы взяли эти точки, не все гармоники периодических возмущений от других планет укладываются четное число раз.

 

 

Рис. 4. Статистическая обработка на форме 6 данных с формы 2 по перигелию Меркурия. Скриншот программы Solsys 7.

 

Так на какое же количество точек в группе нам разбить всю выборку данных. Если мы возьмем вариант разбиения выборки на группы, при котором дисперсия оказалась самая минимальная (критерий min) при 1046 точках в группе, то получим слишком оптимистический результат (572,23+/-0,0009), т.к. возможно, что он получился от случайного стечения обстоятельств. Самым надежным был бы результат, когда и смещение параметра и доверительный интервал рассчитаны по критерию sr (572,16+/-0,79), т.е. как средние значения из всех возможных значений при разных вариантах разбивки первичной выборки на группы. Вот только доверительный интервал по этому критерию получается слишком пессимистический и он значительно больше даже доверительного интервала (+/-0,08), который у нас получился бы, если бы мы статистически обработали все возможные значения полученных смещений параметра в группах. По этому, и этот критерий следует признать не удачным. Тем более что при таком расчете средних значений доверительного интервала не все чисто в статистическом плане.

 

 

Наиболее приемлемым в этом случае критерием для определения значений вековых смещений и предельных отклонений для них я считаю предложенный мною критерий maxmin, т.е. максиминный или минимаксный (в зависимости от того со стороны какого игрока мы смотрим на игру). А говорю я про игру по тому, что эти критерии применяются в теории игр для выбора наилучшей стратегии действий игрока (можете посмотреть здесь [15]). И наилучшей стратегией будет не та, где возможен максимальный выигрыш игрока (в нашем случае минимальная дисперсия), а та где возможен минимальный проигрыш игрока. Вот мы и будем определять максимальный доверительный интервал из минимальных доверительных интервалов так, чтобы все полученные нами доверительные интервалы при разном количестве точек в группах или попадали в него или хотя бы пересекались с ним. А затем уже по этому интервалу найдем матожидание смещения параметра.

 

 

Такой критерий позволяет не только освободить полученные данные по вековым смещениям от дисперсии, вызванной периодическими изменениями параметров орбит, но при этом и гарантирует, что полученный результат не является случайным стечением обстоятельств. Определение вековых смещений параметров и их предельных отклонений с применением этих критериев будет второй частью предложенной мною методики обработки расчетных данных с формы 2. И применять ее мы будем тогда, когда у нас на форме 6 отмечен чекбокс <с перебором всех вариантов разбивки на группы>, по этому, выбирать в рамочке <критерий при переборе> один из рассмотренных выше критериев надо только в том случае, если отмечен этот чекбокс. В противном случае при нажатие на кнопку //Аппроксимировать данные выборки// будет происходить обычная статистическая обработка данных, но с разбивкой всей выборки на группы данных по количеству заданных Вами точек в группе, т.е. будет применяться только первая часть предложенной мною методики обработки таких данных. Рассмотрим определение максиминного значения смещения параметра и его доверительного интервала на примере обработки данных по перигелию Марса (см. рис. 5) по тому, что при обработке данных по Меркурию была слишком сложная для примера схема определения максиминного доверительного интервала.

 

 

 

Рис. 5. Статистическая обработка на форме 6 данных с формы 2 по перигелию Марса. Скриншот программы Solsys 7.

 

Определение верхней и нижней границы максиминного доверительного интервала начинается с того, что за исходный мы принимаем минимальный доверительный интервал из всех возможных. У нас такой интервал получился при разбивке первичных данных на группы по 129 значений и его границы 1599,29 и 1599,95. Теперь проверяем выполняются ли для этого интервала требования предъявляемые нами для максиминного интервала начиная с данных полученных при разбивке первичной выборки по 1 значению, по 2, по 3 и т.д. На рис. 5 мы видим, что это условие выполняется до разбиения выборки по 87 значений, где нижняя граница доверительного интервала (1600,14) получилась больше заданной нами первоначально верхней границы максиминного интервала (1599,95). По этому поднимаем верхнюю границу максиминного интервала до величины 1600,14 и продолжаем сравнивать следующие варианты разбивки первичной выборки на группы уже с новыми границами максиминного интервала. Требования максиминности теперь нарушаются при разбивке выборки по 115 значений в группе. Теперь опять нижняя граница получилась выше верхней границы максиминного интервала. Поднимаем ее до 1600,75 и продолжаем проверку далее. Требование максиминности кругом выполняется и следовательно окончательно нижняя и верхняя границы максиминного интервала будут 1599,29 и 1600,75. Отсюда находим матожидание смещения (1599,29 + 1600,75) / 2 = 1600,02, а доверительный интервал будет  1600,75 - 1600,02 = 0,73, что мы и видим на рис. 5 у надписи максимин.

 

Теперь давайте протестируем предлагаемую мною методику обработки данных по мгновенным значениям параметров орбит получающихся на форме 2 по имитаторам JPL0 и JPL2 (см. табл.6). При этом параметры имитатора JPL0 получены в Лаборатории реактивного движения (JPL, NASA) по данным полученным ими при моделирование (или уже по данным имитатора JPL2), но какова была методика, которую применяли сотрудники JPL для нахождения параметров орбит для имитатора JPL0, мне не известно. Что касается нахождения мною параметров JPL0 maxmin (табл.6), то здесь по всем критериям у меня получаются отличные результаты, т.к. по имитатору JPL0 на форме 2 рассчитывались координаты невозмущенных орбит планет и, следовательно, результаты по определению мною параметров орбит получались очень стабильные и ошибки связаны были только с шагом решения дифференциальных уравнений. Основной шаг решения был 3600 сек, а вблизи перигелия, узла восхождения и при переходе долготы через ноль шаг решения уменьшался. Для Меркурия в 1000 раз, для Венеры в 500, для Земли в 200 и для Марса в 100 раз. При этом данные определялись за период с 1601 года по 2000 год, а для Земли до 1980 года, т.к. вблизи эклиптики из-за периодических возмущений получается очень большой разброс в определение положения узла восхождения, вследствие того, что малейшее отклонение по долготе приводит к резкому смещению точки пересечения орбитой планеты эклиптики, когда две плоскости пересекаются под очень маленьким углом (см. рис. 6). По этому данные по восходящему узлу Земли приходиться обрабатвать только примерно до 1800 года (получается чем меньше этот год, тем больше значение смещения), а чтобы эти данные можно было в нормальном масштабе вывести на график моделирование идет до 1980 года, когда разброс данных не такой большой.

 

 

Таблица 6. Вековые смещения параметров орбит, полученные на программе Solsys7 по различным критериям при переборе данных по всем возможным вариантам разбивки первичной выборки на группы (кроме //JPL2 без перебора// и //JPL0//). Параметры JPL0 maxmin получены по данным с использованием на форме 2 имитатора Солнечной системы JPL0, а остальные имитатора JPL2. (dAlfaP – смещение перигелия, dAlfaU – смещение узла восхождения, dBetta – смещение наклона орбиты, dEks – смещение эксцентриситета, dWsr - изменение угловой скорости и dR - смещение большой полуоси орбиты, а последняя цифра соответствует порядковому номеру планеты).

 

параметр

JPL0

JPL0 maxmin

JPL2 maxmin

JPL2 sr

JPL2 без перебора

dAlfaP1

dAlfaU1

dBetta1

dEks1

dWsr1

dR1

+577,73

-451,22

-21,41

+19,06

0

+55,35

+577,72+/-0,00

-451,23+/-0,00

-21,41+/-0,00

+19,06+/-0,00

0+/-0

+55,0+/-0

+572,20+/-0,07

-449,95+/-0,42

-21,44+/-0,00

+20,50+/-0,06

0+/-0

0+/-0

+572,16+/-0,80

-450,17+/-0,75

-21,44+/-0,02

+20,48+/-0,75

0+/-0

0+/-0

+572,12+/-14,27

-450,27+/-2,56

-21,44+/-0,32

+20,46+/-14,02

0+/-0,02

0+/-0

dAlfaP2

dAlfaU2

dBetta2

dEks2

dWsr2

dR2

+9,66

-999,68

-2,84

-41,07

0

+583,43

+9,65+/-0,00

-999,69+/-0,00

-2,84+/-0,00

-41,07+/-0,00

0+/-0

583,0+/-0

+43,64+/-4,11

-998,36+/-0,33

-2,52+/-0,22

-48,43+/-0,10

0+/-0

0+/-0

+39,32+/-65,30

-998,64+/-1,23

-2,58+/-0,29

-48,39+/-2,03

0+/-0

0+/-0

+40,09+/-1291,05

-998,76+/-10,45

-2,62+/-0,60

-48,39+/-41,69

0+/-0,02

0+/-0

dAlfaP3

dAlfaU3

dBetta3

dEks3

dWsr3

dR3

+1163,77

-

-46,61

-43,92

0

+840,74

+1163,78+/-0,00

-

-46,61+/-0,00

-43,92+/-0,00

0+/-0

840,0+/-0

+1158,55+/-3,34

-861,63+/-6,63

-47,17+/-0,04

-41,93+/-0,17

0+/-0

0+/-0

+1157,75+/-36,59

-849,05+/-77,11

-47,15+/-0,09

-41,87+/-2,83

0+/-0

0+/-0

+1158,29+/-542,29

-831,90+/-694,36

-47,14+/-0,46

-41,99+/-33,47

0+/-0,01

0+/-0

dAlfaP4

dAlfaU4

dBetta4

dEks4

dWsr4

dR4

+1599,88

-1053,25

-29,27

+78,82

0

2763,07

+1599,88+/-0,00

-1053,26+/-0,03

-29,27+/-0,00

+78,82+/-0,00

0+/-0

2763,0+/-0

+1600,02+/-0,73

-1050,49+/-2,36

-28,93+/-0,13

+91,63+/-0,97

0+/-0

0+/-0

+1601,03+/-19,13

-1051,87+/-5,58

-28,99+/-0,22

+91,74+/-8,92

0+/-0

0+/-0

+1601,07+/-299,79

-1052,45+/-22,80

-29,01+/-0,67

+91,17+/-102,05

0+/-0,02

0+/-0,02

 

 

 

Рис. 6. Мгновенные значения угла восходящего узла Земли определенные во время моделирования движения планет Солнечной системы на форме 2. Скриншот программы Solsys7.

 

 

Не значительные различия по смещениям параметров орбит у нас получаются и при обработке данных полученных на имитаторе JPL2, как с перебором всех возможных вариантов разбивки и расчета смещений и доверительных интервалов по разным критериям (maxmin и sr), так и без перебора, т.е. при обычной обработке единичных данных, но вот доверительные интервалы полученных смещений получаются при этом очень разные. По этому, при поиске скорости гравитации я буду использовать при обработке расчетных данных критерий максимина. При этом, в 7-ой версии программы Solsys я уменьшил выборку данных для обработки на форме статистики на две точки, т.е. теперь в отличие от 6-ой версии запись параметров в файл на 2-ой форме идет не со 2-го оборота планеты, а с 3-его (за один оборот параметры определяются два раза). Связано это с тем, что не всегда первоначально заданные параметры орбиты, которые определяют интервал поиска перигелия и узла восхождения, совпадают с получающимися на модели и в результате поиск этих параметров может идти без уменьшения шага решения в нужном интервале углов, что увеличивает погрешность их определения. По этому, если Вы заметите, что приводимые мною сейчас данные немного отличаются от тех, что я рассчитывал по Solsys6 и приводил где либо ранеее, то это объясняется именно этим.

 

 

Теперь что касается интервала времени, который надо задать для проведения вычислительных экспериментов на математической модели. Логично было бы взять интервал с 1800 по 2000 годы, т.е. равный интервалу обработки наблюдательных данных, но при такой маленькой выборке данных у нас получаются значительные доверительные интервалы по этому желательно увеличить этот интервал или использовать другую методику обработки расчетных данных. В таблице 7 представлены данные, обработанные по моей методике по критерию максимина, которые получены на форме 2 программы Solsys7, где моделирование с использованием законов Ньютона производилось с 1401 года по 2000 (у Земли по 1980) с начальными даными рассчитанными по параметрам заданными по теории Ser0 и с шагом решения 3600 секунд, который при определение параметров уменьшался у Меркурия в 1000 раз, у Венеры в 500, у Земли в 200 и у Марса в100. Естественно, при разной продолжительности периода обработки данных и других шагах решения дифференциальных уравнений получатся результаты немного отличные друг от друга, как по смещениям, так и по доверительным интервалам, хотя эти различия будут и не принципиальные.

 

Вот только по восходящему узлу Земли будут очень разные данные, если моделирование начинать с разных моментов времени, когда начальные данные для моделирования рассчитываются с параметрами орбит рассчитанными на эту дату по какой то теории планет. И, если мы выполним моделирование движения объектов Солнечной системы с 1601 года, а не с 1401 года, как дано в таблице 7, то получим для dAlfaU3 -845,52 угл.сек при обработке данных до 1800 года и -818,39 угл.сек до 1900 года, т.е. примерно на 25 угл.сек больше. А, если начать моделирование в 1801 году, то получим всего-навсего -323,40 угл.сек при обработке данных до 1965 года (см. таблицу 9). Ну, а если еще посчитать смещение узла восхождения по воздействию отдельных планет (см. таблицу 8-2), то мы вообще получим положительное смещение +239,90 угл.сек. Таким образом, вследствие не стабильности результатов, этот параметр никак нельзя использовать при расчете критерия оптимизации (также как и перигелий Венеры) при поиске скорости гравитации и мы это и не будем делать задав вес этому параметру в формуле (4-2) равный нулю.

 

 

Таблица 7. Вековые смещения параметров орбит, полученные на программе Solsys7 по критерию  maxmin при обработке расчетных данных за разные периоды времени при начале моделирования в 12 часов 1 января 1401 года.

 

параметр

1801...2000

1601...2000

1401...2000

1401...1800

dAlfaP1

dAlfaU1

dBetta1

dEks1

+529,42+/-0,26

-450,99+/-0,02

-21,44+/-0,01

+20,57+/-0,18

+529,53+/-0,05

-450,20+/-0,43

-21,45+/-0,00

+20,41+/-0,05

+529,62+/-0,01

-449,48+/-1,04

-21,45+/-0,00

+20,47+/-0,03

+529,71+/-0,09

-449,03+/-0,42

-21,45+/-0,00

+20,49+/-0,08

dAlfaP2

dAlfaU2

dBetta2

dEks2

+23,52+/-20,34

-1000,12+/-0,25

-2,77+/-0,00

-49,02+/-0,53

+31,01+/-3,30

-998,91+/-0,11

-2,48+/-0,21

-49,32+/-0,16

+41,98+/-2,22

-997,73+/-1,59

-2,19+/-0,40

-49,46+/-0,04

+48,27+/-1,71

-997,04+/-0,06

-2,03+/-0,18

-49,57+/-0,19

dAlfaP3

dAlfaU3 до 1900

dAlfaU3 до 1800

dBetta3

dEks3

+1146,32+/-9,01

-603,54+/-90,34

-

-47,12+/-0,01

-42,83+/-0,85

+1140,27+/-2,71

-790,29+/-23,42   -821,66+/-14,45

-47,23+/-0,05

-43,29+/-0,10

+1139,16+/-1,70

-823,80+/-3,30

-829,01+/-1,92

-47,30+/-0,09

-43,00+/-0,08

+1139,70+/-2,26

-829,01+/-1,92

-829,01+/-1,92

-47,33+/-0,06

-42,72+/-0,20

dAlfaP4

dAlfaU4

dBetta4

dEks4

+1596,13+/-8,44

-1058,31+/-0,21

-29,20+/-0,02

+96,12+/-1,95

+1596,65+/-1,34

-1051,52+/-3,31

-28,96+/-0,13

+96,87+/-0,47

+1599,17+/-0,44

-1046,74+/-8,16

-28,75+/-0,30

+96,18+/-0,15

+1599,55+/-0,83

-1042,98+/-2,62

-28,63+/-0,11

+96,58+/-0,75

 

 

Как видим с увеличением интервала времени, за который мы обрабатываем данные, прослеживается общая зависимость уменьшения доверительного интервала, но по некоторым смещениям параметров орбит наблюдается наоборот увеличение доверительного интервала, например, для dAlfaU4 и dBetta4. А для некоторых смещений сначала идет уменьшение, а потом увеличение, например, для dAlfaU2. Это объясняется тем, что вековые смещения этих параметров изменяются со временем, т.е. изменение этих параметров со временем имеет не линейный характер. Это можно посмотреть сравнив результаты за периоды с 1401 по 1800 год и с 1601 по 2000 год. По этому, увеличивать интервал времени определения расчетных значений смещений надо до определенного предела, т.е. пока ошибка от нелинейности не станет больше чем от периодических возмущений параметров орбит, но не более чем до 400 лет, т.к. интерполироавать далее 1601 года полученные мною наблюдательные данные для теории Ser0 мы не можем. Но прежде чем переходить к увеличению временного интервала давайте в надежде получить меньшие доверительные интервалы попробуем определить расчетные смещения на интервале в 200 лет как сумму смещений при возмущение нашей исследуемой планеты только одной другой планетой, как мы это делали для перигелия Меркурия в табл. 1.

 

 

Таблица 8-1. Расчетные значения смещения перигелиев орбит планет (в угл.сек. за 100 лет) от действия на них отдельных планет с системой масс JPL за период времени с 1801 по 2000 год.

 

 

Меркурий

Венера

Земля

Марс

Меркурий

+0,006+/-0,003

-147,945+/-0,204

-13,585+/-0,099

+0,763+/-0,028

Венера

+275,950+/-0,091

+0,002+/-0,003

+341,837+/-5,898

+49,881+/-0,942

Земля

+90,196+/-0,023

-558,944+/-11,995

+0,003+/-0,006

+228,374+/-1,460

Марс

+2,467+/-0,003

+74,696+/-6,729

+99,476+/-3,495

+0,002+/-0,007

Юпитер

+153,545+/-0,174

+674,666+/-11,611

+705,877+/-15,863

+1247,544+/-4,50

Сатурн

+7,322+/-0,058

+8,184+/-0,936

+19,361+/-1,590

+66,301+/-0,863

Уран

+0,152+/-0,009

+0,236+/-0,061

+0,562+/-0,019

+1,082+/-0,081

Нептун

+0,036+/-0,004

+0,047+/-0,175

+0,078+/-0,062

+0,493+/-0,030

Плутон

+0,006+/-0,002

+0,002+/-0,002

+0,003+/-0,006

+0,002+/-0,007

итого

+529,680+/-0,376

+50,944+/-31,716

+1153,612+/-27,04

+1594,442+/-7,92

 

 

Таблица 8-2. Расчетные значения смещения узлов восхождения планет (в угл.сек. за 100 лет) от действия на них отдельных планет с системой масс JPL за период времени с 1801 по 2000 год.

 

 

Меркурий

Венера

Земля*

Марс

Меркурий

+0,002+/-0,002

+11,227+/-0,005

+615,311+/-0,621

+1,819+/-0,010

Венера

-194,064+/-0,011

+0,001+/-0,003

+10805,54+/-224,68

+30,339+/-0,194

Земля

-99,551+/-0,010

-724,176+/-0,164

-0,009+/-0,008

-225,531+/-0,264

Марс

-1,924+/-0,004

-4,717+/-0,018

+1207,999+/-1,088

-0,001+/-0,008

Юпитер

-148,196+/-0,037

-273,106+/-0,073

-10805,104+/-74,01

-833,697+/-0,312

Сатурн

-6,984+/-0,012

-8,301+/-0,011

-1578,547+/-2,510

-27,080+/-0,113

Уран

-0,136+/-0,003

-0,276+/-0,006

+3,033+/-0,144

-0,754+/-0,045

Нептун

-0,048+/-0,005

-0,008+/-0,004

-8,303+/-0,019

-0,330+/-0,018

Плутон

+0,002+/-0,002

+0,001+/-0,003

-0,015+/-0,004

-0,001+/-0,008

итого

-450,899+/-0,086

-999,431+/-0,287

+239,903+/-303,139

-1055,236+/-0,97

 

* - данные для Земли от Венеры и Юпитера получены на интервале от 1801 до 1870 года, т.к. дальше смещение от них носит явно выраженный не линейный характер, но и при этом ограничении результат получился не удовлетворительным (получилось даже положительное значение), т.к. влияние этих двух планет на смещение восходящего узла Земли просто огромно по сравнению с другими планетами.

 

 

 

 

 

 

Таблица 8-3. Расчетные значения смещения углов наклона орбит планет (в угл.сек. за 100 лет) от действия на них отдельных планет с системой масс JPL за период времени с 1801 по 2000 год.

 

 

Меркурий

Венера

Земля

Марс

Меркурий

+0+/-0

+1,183+/-0,001

-0,276+/-0,000

+0,009+/-0,000

Венера

-14,788+/-0,003

+0+/-0

-28,887+/-0,057

-1,290+/-0,005

Земля

-1,360+/-0,004

+0,250+/-0,009

+0+/-0

+0,112+/-0,001

Марс

-0,029+/-0,000

+0,133+/-0,001

-0,788+/-0,001

0+/-0

Юпитер

-4,890+/-0,009

-3,759+/-0,005

-15,565+/-0,026

-25,520+/-0,011

Сатурн

-0,420+/-0,002

-0,522+/-0,001

-1,240+/-0,001

-2,466+/-0,003

Уран

-0,003+/-0,000

+0+/-0

-0,008+/-0,000

-0,006+/-0,000

Нептун

-0,002+/-0,000

-0,003+/-0,000

-0,004+/-0,000

-0,010+/-0,000

Плутон

+0+/-0

+0+/-0

+0+/-0

0+/-0

итого

-21,492+/-0,017

-2,716+/-0,017

-46,788+/-0,085

-29,171+/-0,020

 

 

Таблица 8-4. Расчетные значения смещения эксцентриситетов орбит планет (увеличино в 1 000 000 раз и за 100 лет) от действия на них отдельных планет с системой масс JPL за период времени с 1801 по 2000 год.

 

 

Меркурий

Венера

Земля

Марс

Меркурий

+0+/-0

-7,843+/-0,007

-0,701+/-0,003

+0,203+/-0,002

Венера

+13,483+/-0,100

+0+/-0

+7,011+/-0,266

+0,348+/-0,350

Земля

+5,511+/-0,020

-24,080+/-0,474

+0+/-0

+10,275+/-0,725

Марс

-0,294+/-0,000

-1,107+/-0,223

-7,499+/-0,543

0+/-0

Юпитер

+1,582+/-0,239

-15,343+/-0,428

-40,702+/-0,958

+77,324+/-1,905

Сатурн

+0,285+/-0,051

-0,338+/-0,067

-0,268+/-0,171

+3,140+/-0,392

Уран

+0,005+/-0,005

-0,004+/-0,001

-0,006+/-0,001

+0,061+/-0,021

Нептун

-0,010+/-0,011

-0,009+/-0,002

-0,019+/-0,011

+0,008+/-0,021

Плутон

+0+/-0

+0+/-0

+0+/-0

0+/-0

итого

+20,561+/-0,424

-48,720+/-1,202

-42,184+/-1,953

+91,359+/-3,416

 

 

Как видим результаты получились гораздо хуже, чем при определение доверительных интервалов, когда на исследуемую планету действовали сразу все возмущения от остальных планет. Таким образом увеличить точность определения параметров орбит по предлагаемой мною методике обработки данных, получающихся при моделировании движения планет, можно только соответствующим подбором интервала времени, за который будут определяться смещения параметров орбит. А чтобы при поиске скорости гравитации результаты были сопоставимы, я все эксперименты буду делать при одинаковой продолжительности экспериментов и при тех же шагах решения. Сейчас я считаю, что наиболее приемлемым является промежуток времени от 200 до 400 лет, т.е. с 1601 или 1801 по 2000 год. Но, т.к. наблюдаемые смещения в теории Ser0 были получены мною при обработке данных примерно с 1800 по 1980 год, то и при поиске скорости гравитации желательно проводить вычислительные эксперименты именно в этом интервале лет, не смотря на то, что доверительные интервалы полученных смещений параметров орбит при этом будут хуже, чем при интервале с 1600 по 2000 год. По этому я сейчас приведу таблицу 9, где не только дам получившиеся у меня наблюдаемые смещения параметров орбит для теории Ser0, но и получающиеся при этом по теории Ньютона их расчетные значения на интервале с 1801 по 2000 год. А, чтобы еще раз показать преимущества моей методики обработки данных по критерию maxmin над стандартной, приведу данные, полученные при обработке данных полученных при моделировании, которые дает как одна так и вторая методики.

 

 

Таблица 9. Вековые смещения параметров орбит и их доверительные интервалы, полученные на программе Solsys7 с использованием математической модели Солнечной системы, в которой используются законы Ньютона, (в интервале с 1801 по 2000 год), когда статистическая обработка данных ведется по критерию maxmin с перебором различных вариантов разбивки данных на группы, и классически по одной точке, а также полученные на программе Solsys7 при обработке непосредственно данных оптических наблюдений для создания кинематической астрономической теории Ser0 (в интервале примерно с 1800 по 1980 год) при классической статистической обработке данных.

 

параметры

теория Ser0

модель

maxmin

классика

dAlfaP1 угл.сек. /100 лет

dAlfaU1 угл.сек. /100 лет

dBetta1 угл.сек. /100 лет

dEks1 * 1 000 000

+578,04 +/- 7,25

-433,15 +/ 8,21

-19,84 +/- 0,51

+20,10 +/- 3,98

+529,71 +/- 0,27

-451,41 +/ 0,01

-21,45 +/- 0,01

+20,49 +/- 0,20

+529,68 +/- 14,82

-451,53 +/- 2,38

-21,45 +/- 0,34

+20,46 +/- 13,42

dAlfaP2

dAlfaU2

dBetta2

dEks2

+30,37 +/- 83,48

-996,50 +/- 6,84

-2,83 +/- 0,21

-45,32 +/- 2,85

+57,73 +/- 24,91

-1000,64 +/- 0,13

-2,82 +/- 0,00

-48,87 +/- 0,71

+59,00 +/- 1307,90

-1000,81 + 10,60

-2,87 +/- 0,52

-48,78 +/- 41,57

dAlfaP3

dAlfaU3

dBetta3

dEks3

+1141,58 +/- 9,24

-869,00 +/- 23,90

-47,45 +/- 0,13

-41,57 +/- 0,82

+1159,43 +/- 2,53

-323,40 +/- 26,78

-47,16 +/- 0,01

-43,17 +/- 0,75

+1162,88 +/- 560,09

-18,24 +/- 1976,78

-47,09 +/- 0,47

-42,60 +/- 33,00

dAlfaP4

dAlfaU4

dBetta4

dEks4

+1582,74 +/- 30,23

-1028,02 +/- 29,93

-29,39 +/- 0,85

+95,53 +/- 12,37

+1588,01 +/- 9,52

-1060,25 +/- 0,54

-29,31 +/- 0,01

+93,07 +/- 2,22

+1593,18 +/- 301,46

-1060,60 +/- 22,52

-29,33 +/- 0,64

+90,16 +/- 102,25

 

 

Но относительно большие доверительные интервалы по некоторым параметрам даже при обработке данных по критерию maxmin это еще не все неприятности, которые будут мешать нам при поиске скорости гравитации. Оказывается, на результат будет влиять даже теория планет, по которой мы будем задавать начальные данные для начала моделирования на заданный момент времени. И в таблице 10 я привожу получившиеся у меня результаты при задание начальных данных по теории JPL2  и по теории Ser0. При этом я привожу для расчета параметров орбит при одинаковых начальных данных (рассчитанных по теории Ser0) два варианта моделирования - когда используется закон притяжения Ньютона и закон Ньютона с поправкой Холла (при скорости гравитации равной бесконечности). И получается, что полученные данные по смещениям параметров орбит (кроме перигелиев) при расчете начальных данных по теории Ser0 даже больше отличаются от полученных данных по смещениям параметров орбит, когда начальные данные задавали по теории JPL2, а движение планет в обоих случаях рассчитывали по одному и тому же закону Ньютона, чем от полученных смещений параметров орбит когда задавали начальные данные по той же теории - Ser0, но при расчете движения планет уже по теории Холла.

 

Таблица 10. Вековые смещения параметров орбит, полученные на программе Solsys7 по критерию  maxmin при обработке рассчетных данных за период времени с 1601 по 2000 годы при разных начальных данных для моделирования и с использованием системы масс JPL для JPL2 и Ser0. Для начальных данных рассчитанных по теории Ser0 выполнено два варианта расчета – по закону Ньютона и по этому закону с поправкой Холла.

 

параметр

JPL2 – R^2

Ser0 – R^2

Ser0 – R^n*

dAlfaP1

dAlfaU1

dBetta1

dEks1

+529,67+/-0,09

-450,39+/-0,43

-21,456+/-0,003

+20,424+/-0,066

+529,59+/-0,08

-450,36+/-0,42

-21,455+/-0,002

+20,421+/-0,073

+573,41+/-0,07

-450,40+/-0,44

-21,457+/-0,002

+20,417+/-0,073

dAlfaP2

dAlfaU2

dBetta2

dEks2

+44,789+/-4,631

-999,21+/-0,34

-2,530+/-0,222

-49,511+/-0,127

+45,555+/-4,331

-999,08+/-0,359

-2,521+/-0,222

-49,417+/-0,133

+61,193+/-4,379

-999,06+/-0,33

-2,521+/-0,222

-49,409+/-0,132

dAlfaP3

dAlfaU3**

dBetta3

dEks3

+1142,66+/-3,25

-828,54+/-16,85

-47,240+/-0,046

-43,262+/-0,226

+1141,97+/-3,071

-845,52+/-12,34

-47,235+/-0,045

-43,291+/-0,224

+1152,56+/-3,05

-845,36+/-12,43

-47,234+/-0,044

-43,291+/-0,225

dAlfaP4

dAlfaU4

dBetta4

dEks4

+1598,00+/-0,52

-1052,24+/-2,73

-29,009+/-0,131

+96,245+/-0,708

+1598,09+/-0,66

-1052,79+/-2,72

-29,025+/-0,132

+96,255+/-0,714

+1603,70+/-0,60

-1052,77+/-2,76

-29,025+/-0,132

+96,252+/-0,713

 

*  - показатель степени n принят таким, каким его брал Ньюком, т.е. n=2,0000001612, а не 2,0000001574, как рассчитал его Холл.

** - dAlfaU3 рассчитывалось за интервал с 1601 по 1800 год.

 

 

Теперь непосредственно о задание начальных данных для начала процесса моделирования. Что касается теории JPL2, т.е. эфемерид DE405, то здесь никаких проблем, т.к. по ним мы как раз и вычисляем мгновенные значения координат и скоростей планет по осям X, Y и Z, т.е. находим то, что получилось у сотрудников JPL на их модели. Причем координаты и скорости будут при этом заданы относительно барицентра Солнечной системы, который мы можем просто поместить в начало системы координат. А вот, когда мы рассчитываем начальные данные по какой то аналитической теории планет, то здесь у нас возникают проблемы. В теориях планет дается расчет параметров орбит планет при их движение относительно неподвижного Солнца (аналогично и для Луны) и координаты планет мы можем рассчитать по этим параметрам по известным формулам, которые приведены в любом учебнике по астрономии. Расписывать эти формулы я не буду, а просто приведу код подпрограммы для их расчета в программе Solsys7.

 

Public Sub RasKoordinat() 'расчет координат без периодических возмущений

AlfaM0(i) = FRAC2(AlfaL0(i) / 2 / pi) * 2 * pi - AlfaP0(i)

dUgolE = 1 * 10 ^ -10 '                   производить итерации до этой точности

If AlfaM0(i) < 0 Then AlfaM0(i) = AlfaM0(i) + 2 * pi

UgolE0(0) = AlfaM0(i)

    For k = 1 To 99

    UgolE0(k) = AlfaM0(i) + Eks0(i) * Sin(UgolE0(k - 1))

    If Abs(UgolE0(k - 1) - UgolE0(k)) > dUgolE Then GoTo 100

    UgolE = UgolE0(k): k = 100 '                эксцентрическая аномалия

100: Next k

 

UgolV(i) = 2 * Atn(Sqr((1 + Eks0(i)) / (1 - Eks0(i))) * Tan(UgolE / 2)) 'истинная аномалия

Rsr(i) = R0sr(i) * (1 - Eks0(i) * Eks0(i)) / (1 + Eks0(i) * Cos(UgolV(i)))

UgolU(i) = UgolV(i) + AlfaP0(i) - AlfaU0(i) '           аргумент широты

X(i) = Rsr(i) * (Cos(UgolU(i)) * Cos(AlfaU0(i)) - Sin(UgolU(i)) * Sin(AlfaU0(i)) * Cos(Betta0(i)))

Y(i) = Rsr(i) * (Cos(UgolU(i)) * Sin(AlfaU0(i)) + Sin(UgolU(i)) * Cos(AlfaU0(i)) * Cos(Betta0(i)))

Z(i) = Rsr(i) * Sin(UgolU(i)) * Sin(Betta0(i))

'для Луны это будут координаты от Земли, а для планет от Солнца

End Sub

 

 

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

 

Public Sub RasSkorostey() 'расчет скоростей без периодических возмущений

Ipl0 = 0: If kodMoon = 1 Then Ipl0 = 13 '13 - это Земля без Луны

m(12) = m(Ipl0): If kodBary = 1 Then m(12) = m(Ipl0) + m(i)

'kodBary = 0 - расчет относительно центрального тела

'kodBary = 1 - расчет относительно барицентра двух тел

'kodMoon = 0 - если определяем скорость планеты около Солнца

'kodMoon = 1 - если определяем скорость Луны около Земли

p = R0sr(i) * (1 - Eks0(i) * Eks0(i))

Vr = Sqr(gamma * m(12) / p) * Eks0(i) * Sin(UgolV(i))

Vn = Sqr(gamma * m(12) / p) * (1 + Eks0(i) * Cos(UgolV(i)))

VX(i) = X(i) * Vr / Rsr(i) + Vn * (-Sin(UgolU(i)) * Cos(AlfaU0(i)) - Cos(UgolU(i)) * Sin(AlfaU0(i)) * Cos(Betta0(i)))

VY(i) = Y(i) * Vr / Rsr(i) + Vn * (-Sin(UgolU(i)) * Sin(AlfaU0(i)) + Cos(UgolU(i)) * Cos(AlfaU0(i)) * Cos(Betta0(i)))

VZ(i) = Z(i) * Vr / Rsr(i) + Vn * Cos(UgolU(i)) * Sin(Betta0(i))

End Sub

 

 

К счастью для планет Солнечной системы эта ошибка от некорректного расчета начальных скоростей будет не очень заметна, а вот для Луны, масса которой всего в 81 раз меньше массы Земли, это будет уже очень заметно. По этому для рассмотрения вопроса различных вариантов расчета начальных данных я создал даже дополнительную форму 20 (Расчет начальных данных для Луны), скриншот которой Вы видите на рис. 7. Полностью все возможности, которые предоставляет эта форма, Вы прочитаете в описание программы, а сейчас только замечу, что на ней Вы можете подкорректировать полученные при стандартном расчете, как координаты, так и скорости Земли и Луны, чтобы минимизировать ошибку от применения стандартных формул расчета. Например, на рисунке мы видим расчет начальных координат и скоростей Земли и Луны, когда Земля находится в центре эллипса, по которому Луна вращается вокруг Земли. При таком стандартном расчете мы имеем ошибку по перигелию в -7,1 градуса и по большой полуоси в -4,2 тыс.км. (остальные параметры орбиты определенные при моделирование совпали с начальными или точно в ноль или с минимальной ошибкой). Но, если сделать предлагаемые мною уточнения координат и скоростей Земли и Луны, то, во-первых, не будет движения барицентра, как мы это видим на рисунке, и пропадет ошибка по перигелию, но по большой полуоси все-таки останется +4,7 тыс.км и в добавок теперь не будет сохраняться количество движения системы по оси Х, хотя энергия системы и момент количества движения системы будут продолжать сохраняться. А, если сделать стандартный расчет для вращения Луны вокруг барицентра, то все ошибки по получающимся параметрам орбит будут совпадать с начальными или точно в ноль или с минимальной погрешностью, но барицентр будет двигаться, как и при первом стандартном расчете. Но, если сделать и в этом расчете предлагаемые мною уточнения, то движение барицентра прекратится и все погрешности в параметрах будут минимальные или нулевые при неизменных энергии, моменте количества движения и количестве движения по оси Х.

 

 

Рис. 7. Варианты расчета начальных данных для начала моделирования для Луны и Земли по параметрам орбиты Луны на форме 20. Скриншот программы Solsys7.

 

Вот только не все так прекрасно, т.к., во-первых, мои поправки не все строго теоретически обоснованы, а, во-вторых, второй вариант расчета имеет смысл только для одной планеты, а не для всей Солнечной системы. По этому, когда будете задавать на 1-ой форме параметры орбит и рассчитывать по ним начальные данные для моделирования, то можете, конечно, отметить чекбокс //Bary//, чтобы расчет велся по второму варианту, но учтите, что поправки, приведенные на форме 20 будут при этом будут учитываться только для Луны и Земли, если Луна при моделировании будет рассматриваться отдельно от Земли. Ведь при таком расчете координат и скоростей планеты относительно барицентра Вы будете учитывать только координаты барицентра этой планеты и Солнца, а на самом деле у всей Солнечной системы барицентр будет в другом месте, т.е. Солнце будет находиться в другом месте. В таком варианте общий расчет для всех планет очень сильно усложнится, по этому, этот усложненный вариант расчета для планет я даже не рассматривал. Ведь с такими же погрешностями можно рассчитать координаты и скорости планет и при движение относительно Солнца (не отмечать чекбокс //Bary//). Потом можете для обоих вариантов найти скорость Солнца и на нее уточнить скорости планет и найти координаты барицентра системы и уточнить координаты планет. В любом случае, как я писал, для планет это не критично и можете считать и так и так. Даже при расчете относительно Солнца получится, наверное, точнее.

 

 

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

 

 

А теперь давайте перейдем к аномальным остаткам смещения параметров орбит, а для того, чтобы знать на изменение каких параметров при этом обращать особое внимание, я определю аномальные остатки смещения параметров орбит и их доверительные интервалы. При этом полученные мною в предыдущей статье значения вековых изменений параметров Wsr и R я не включил в табл. 5, т.к. их следует признать не пригодными для поиска скорости распространения гравитации, как не достоверные и по этому для них я не буду находить и аномальные остатки. Но прежде чем перейти к получающимся у меня аномальным остаткам, т.е. остаткам, которые не объясняются теорией Ньютона, давайте посмотрим какие были аномальные остатки у Ньюкома (табл. 11-1). Он приводит аномальные остатки вычисленные в текущей эклиптике, т.е. в эпохе даты, и по этому у него отсутствует понятие наклона орбиты Земли к плоскости эклиптики и не может быть по определению восходящего узла орбиты Земли. Но он все же приводит для Земли изменение угла наклона, но не орбиты Земли к фиксированной эклиптике, а эклиптики к экватору. В своей таблице он также приводит не смещение перигелия, а его произведение на эксцентриситет планеты и не смещение восходящего узла, а его произведение на синус угла наклона орбиты, но нам нужны именно смещения параметров орбит, по этому я в табл. 11-1 привожу данные Ньюкома по перигелию уже деленные на эксцентриситеты планеты и угол восходящего узла деленый на угол наклона. Эти значения эксцентриситетов планет и углов наклона я взял из теории Ньюкома, которые получаются для 1900 года. А его значения смещения эксцентриситета я перевел из угл.сек. в безразмерные величины.

 

 

Да действительно мы видим четыре аномалии, о которых говорил Ньюком (перигелии Меркурия и Марса, эксцентриситет Меркурия и узел Венеры), но мы видим и еще три аномалии – узел Меркурия, угол наклона Венеры и эксцентриситет Марса. А вот почему Ньюком их не посчитал аномальными не понятно. В своей работе [6] он пишет, что у упомянутых им четырех аномалий отклонения превышали свои вероятные ошибки, но в таблице он приводит средние ошибки (как я понимаю среднеквадратичное отклонение, т.е. стандарт). Но для получения вероятной ошибки надо среднюю ошибку умножить на 0,67454, а в таком случае не только смещение узла Меркурия, угла наклона Венеры и эксцентриситет Марса остаются аномальными, но по этому критерию следует и смещение эксцентриситета Венеры и перигелия и угла наклона Земли тоже признать аномальными. В общем, тут есть над чем подумать и в частности над тем как Ньюком определял суммарные ошибки (у него они не равны сумме наблюдаемых и расчетных), а пока я возьму из этой таблицы значение расчетных ошибок, которое буду считать среднеквадратичным отклонением при определение своих значений аномальных остатков, которые буду считать таковыми если отклонение наблюдаемого смещения параметра от расчетного будет превышать сумму доверительных интервалов определения наблюдаемых и расчетных значений.

 

Таблица 11-1. Вековые смещения параметров орбит и их аномальные остатки с указанием доверительного интервала (ошибок) при доверительной вероятности (надежности) 90% по данным Ньюкома [6]. (смещения углов даны в угл.сек., а эксцентриситета увеличены в 1 000 000 раз).

 

 

параметры

Смещения параметров

наблюдаемые

расчетные

остаток

dAlfaP1

dAlfaU1

dBetta1

dEks1

575,06+/-1,95

-753,70+/-3,69

7,14+/-0,80

16,29+/-2,42

533,82+/-0,78

-758,70+/-1,31

6,76+/-0,01

20,56+/-0,05

+41,24+/-2,09

+5,00+/-4,27

+0,38+/-0,80

-4,27+/-2,42

dAlfaP2

dAlfaU2

dBetta2

dEks2

42,52+/-29,33

-1780,56+/-2,03

3,87+/-0,30

-45,86+/-0,97

49,85+/-21,99

-1790,69+/-2,03

3,49+/-0,14

-46,88+/-1,16

-7,33+/-36,66

+10,14+/-2,87

0,38+/-0,33

+1,02+/-1,50

dAlfaP3

dAlfaU3

dBetta3

dEks3

1162,92+/-7,16

-

-47,11+/-0,23

-41,45+/-0,44

1156,95+/-2,98

-

-46,89+/-0,09

-41,55+/-0,19

+5,97+/-7,76

-

-0,22+/-0,27

+0,10+/-0,48

dAlfaP4

dAlfaU4

dBetta4

dEks4

1602,69+/-3,75

-2248,44+/-6,19

-2,26+/-0,20

92,11+/-1,31

1594,65+/-0,43

-2249,37+/-2,79

-2,25+/-0,04

90,71+/-0,05

+8,04+/-3,75

+0,93+/-6,81

-0,01+/-0,20

+1,4+/-1,31

 

 

При проведение научных исследований доверительная вероятность (надежность) данных должна быть 95%, а при такой надежности данных доверительные интервалы должны браться как два стандарта (sigma), т.е. как два среднеквадратичных отклонения. И для наблюдательных данных я приму ошибки их определения такими как получились у меня доверительные интервалы при 95% надежности данных в предыдущей статье [4] при создание теории Ser0, а ошибки Ньюкома для наблюдательных данных в табл. 11-2 и 11-3 я привожу только для справки. При этом обращаю Ваше внимание на то, что в статье [4] я приводил предельные отклонения для самих параметров орбит планет на интервале примерно 200 лет, которые получаются при стандартной статистической обработке данных на форме 6 программы Solsys7, а в таблицах 11-2 и 11-3 я привожу предельные отклонения для смещений параметров орбит. В связи с этим приведенные в работе [4] предельные отклонения параметров орбит я уменьшил в 2 раза, т.к. , если, например, для самого параметра - перигелия Меркурия предельное отклонение было +/-7,25 угл.сек, то для его векового смещения, т.е. за 100 лет оно будет только +/-3,62 угл.сек, а при этом за 200 лет ошибка по самому параметру и получится +/-7,25 угл.сек.

 

 

А вот расчетные ошибки у Ньюкома это ошибки связанные с ошибками определения масс планет, а мои расчетные ошибки появляются из за обработки мгновенных значений параметров орбит определенных при моделировании. У меня данных по ошибкам связанным с разными системами масс планет нет, по этому эти ошибки я приму такими как их определил Ньюком. При этом у Ньюкома не было статистических ошибок при определение расчетных смещений параметров орбит, т.к. он определял их аналитически при воздействие через пертурбационную функцию всех планет на одну исследуемую планету. А предложенная мною методика обработки мгновенных параметров орбит найденных при моделировании для получения вековых смещений по критерию максимина дает значения смещений только в доверительном интервале при 95% надежности и по этому суммарные расчетные ошибки я буду брать как сумму ошибок Ньюкома и своих. При этом, в табл. 11-2 и 11-3 я приму свои расчетные ошибки, которые у меня получились в табл. 9 для периода 1801...2000 годы

Таблица 11-2. Вековые смещения параметров орбит и их аномальные остатки с указанием доверительного интервала (ошибок) при доверительной вероятности (надежности) 95% полученные мною при обработке расчетных и наблюдаемых данных за 200 лет. (смещения углов даны в угл.сек., а эксцентриситета увеличены в 1 000 000 раз)

 

 

параметры

Ошибки определения

Смещения параметров

наблюдения

расчет

New*

Ser

New*

Ser

наблюдения

расчет

остаток

dAlfaP1

dAlfaU1

dBetta1

dEks1

3,90

7,38

1,6

4,84

7,25

8,21

0,51

3,98

1,56

2,62   0,02

0,10

0,27

0,01    0,01

0,20

578,04+/-7,25

-433,15+/-8,21

-19,84+/-0,51

20,10+/-3,98

529,71+/-1,83

-451,41+/-2,63

-21,45+/-0,03

20,49+/-0,30

+48,33+/-9,08

+18,26+/-10,84

+1,61+/-0,54

-0,39+/-4,28

dAlfaP2

dAlfaU2

dBetta2

dEks2

58,66

4,06

0,60

1,94

83,48

6,84

0,21

2,85

43,98

4,06   0,28

2,32

24,91

0,13     0,00

0,71

30,37+/-83,48

-996,50+/-6,84

-2,83+/-0,21

-45,32+/-2,85

57,73+/-78,89

-1000,64+/-4,19

-2,82+/-0,28

-48,87+/-3,03

-27,36+/-162,37

+4,14+/-11,03

-0,01+/-0,49

+3,55+/-5,88

dAlfaP3

dAlfaU3

dBetta3

dEks3

14,32

-

0,46

0,88

9,24

-

0,13

0,82

5,96

-        0,18

0,38

2,53

26,78    0,01

0,75

1141,58+/-9,24

-

-47,45+/-0,13

-41,57+/-0,82

1159,43+/-8,49

-323,40+/-26,78

-47,16+/-0,19

-43,17+/-1,13

-17,85+/-17,73

-

-0,29+/-0,32

+1,60+/-1,95

dAlfaP4

dAlfaU4

dBetta4

dEks4

7,5

12,38

0,40

2,62

30,23

29,93

0,85

12,37

0,86

5,58      0,08

0,10

9,52

0,54    0,01

2,22

1582,74+/-30,23

-1028,02+/-29,93

-29,39+/-0,85

95,53+/-12,37

1588,01+/-10,38

-1060,25+/-6,12

-29,31+/-0,09

93,07+/-2,32

-5,27+/-40,61

+32,23+/-36,05

-0,08+/-0,94

+2,46+/-14,69

 

* - расчетные ошибки Ньюкома увеличены мною в этой таблице в два раза, чтобы доверительная вероятность и его данных была 95% как у меня, а не 90% как было у него при средних ошибках, а ошибки по наблюдаемым данным даны справочно.

 

Как видим, из 7-и аномальных остатков вековых смещений параметров орбит по данным Ньюкома, у меня остались только два – перигелий и узел Меркурия, но к ним добавился еще угла наклона Меркурия. А, если уменьшить доверительную вероятность моих наблюдаемых и расчетных данных до 90%, т.е. сделать ее такой, какой она была у Ньюкома, при определение им аномальных остатков, когда доверительный интервал будет равен одному стандарту, то у нас аномальных остатков будет больше (см. табл.11-3), но это не радует. Ведь по большому счету ярко выраженными аномальными остатками можно считать только три остатка у Меркурия, а для поиска скорости гравитации надо, чтобы были ярко выраженных хотя бы два аномальных остатка, но у разных планет. Но ничего не поделаешь. Имеем то, что имеем. По этому, поиск скорости гравитации будем вести по аномальным остаткам всех параметров всех планет, но при этом будем задавать разные весовые коэффициенты полученным остаткам для их суммарного учета в едином критерии оптимизации при проведении вычислительных экспериментов по планам многофакторного планирования. И чем более достоверен получившийся в табл. 11-3 аномальный остаток тем будет больше его вес в расчете критерия оптимизации. Хотя возможен и такой вариант, что весовые коэффициенты я буду выставлять по процентному отношению ошибки определения наблюдаемых смещений к величине этих смещений. Надеюсь, что суммарное влияние скоростей гравитации и Солнечной системы на критерий оптимизации, рассчитанный по изменению всех смещений, позволит мне все же найти оптимальную скорость гравитации по имеющимся у меня наблюдаемым смещениям параметров орбит.

 

 

 

 

Таблица 11-3. Вековые смещения параметров орбит и их аномальные остатки с указанием доверительного интервала (ошибок) при доверительной вероятности (надежности) 90% полученные мною при обработке расчетных и наблюдаемых данных за 200 лет. (смещения углов даны в угл.сек., а эксцентриситета увеличены в 1 000 000 раз)

 

 

параметры

Ошибки определения

Смещения параметров

наблюдения

расчет

New*

Ser**

New*

Ser

наблюдения

расчет

остаток

dAlfaP1

dAlfaU1

dBetta1

dEks1

1,95

3,69

0,80

2,42

3,63

4,11

0,26

1,99

0,78

1,31   0,01

0,05

0,14

0,00    0,00

0,10

578,04+/-3,63

-433,15+/-4,11

-19,84+/-0,26

20,10+/-1,99

529,71+/-0,92

-451,41+/-1,31

-21,45+/-0,01

+20,49+/-0,15

+48,33+/-4,55

+18,26+/-5,42

+1,61+/-0,27

-0,39+/-2,14

dAlfaP2

dAlfaU2

dBetta2

dEks2

29,33

2,03

0,30

0,97

41,74

3,42

0,11

1,43

21,99

2,03   0,14

1,16

12,46

0,07     0,00

0,36

30,37+/-41,74

-996,50+/-3,42

-2,83+/-0,11

-45,32+/-1,43

57,73+/-34,45

-1000,64+/-2,10

-2,82+/-0,14

-48,87+/-1,52

-27,36+/-76,19

+4,14+/-5,52

-0,01+/-0,25

+3,55+/-2,95

dAlfaP3

dAlfaU3

dBetta3

dEks3

7,16

-

0,23

0,44

4,62

-

0,07

0,41

2,98

-        0,09

0,19

1,27

13,39    0,00

0,38

1141,58+/-4,62

-

-47,45+/-0,07

-41,57+/-0,41

1159,43+/-4,25

-323,40+/-13,39  -47,16+/-0,09

-43,17+/-0,57

-17,85+/-8,87

-

-0,29+/-0,16

+1,60+/-0,98

dAlfaP4

dAlfaU4

dBetta4

dEks4

3,75

6,19

0,20

1,31

15,12

14,97

0,43

6,19

0,43

2,79      0,04

0,05

4,76

0,27    0,00

1,11

1582,74+/-15,12

-1028,02+/-14,97

-29,39+/-0,43

95,53+/-6,19

1588,01+/-5,19

-1060,25+/-3,06

-29,31+/-0,04

93,07+/-1,16

-5,27+/-20,31

+32,23+/-18,03

-0,08+/-0,47

+2,46+/-7,35

 

* - расчетные ошибки Ньюкома даны с доверительной вероятность 90% , т.е. как было у него при средних ошибках, а ошибки по наблюдаемым данным даны справочно.

** - мои ошибки (доверительные интервалы) полученные для наблюдаемых параметров орбит в теории  Ser0 я уменьшил в 2 раза, т.к. там они определялись мною с доверительной вероятность 95% и доверительный интервал брался тогда 2 сигма, а с доверительной вероятность 90%  доверительный интервал будет примерно 1 сигма, хотя это и не строго, т.к. я применял минимаксную методику обработки данных, а не стандартную.

 

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

 

Y = k0 + k1*X1 + k2*X2 + k3*X3 + k4*X4 +

+ k5*X1*X2 + k6*X1*X3 + k7*X1*X4 + k8*X2*X3 + k9*X2*X4 + k10*X3*X4 +

+ k11*X1^2 + k12*X2^2 + k13*X3^2 + k14*X4^2                                     (3)

 

где Y – критерий оптимизации, который надо минимизировать или максимизировать, X1, X2, X3, X4 – оптимизируемые параметры (факторы), а k0...k14 – коэффициенты, которые мы получаем методом наименьших квадратов при статистической обработке экспериментальных данных.

 

Специфическим моментом, связанным с теорией многофакторного планирования, является то, что оптимизацию ведут почти всегда по отклику системы (в нашем случае это будет по самим смещениям параметров орбит планет, а не по отклонению расчетных смещений от наблюдаемых) и находят оптимальные параметры, когда этот отклик (критерий оптимизации) принимает максимальное или минимальное значение. Но нам не надо найти параметры системы, когда наши вековые смещения параметров орбит будут максимальные или минимальные. Нам надо, чтобы аномальный остаток между наблюдаемыми смещениями и рассчитанными при моделировании был минимальным. Да, с практическим применением этого критерия (см. формулу 4-4) Вы уже знакомы по предыдущей статье, когда мы оптимизировали параметры орбит планет. Но тогда я не стал приводить никаких обоснований возможности применения этого критерия, а просто сослался на эту статью, которая еще не была написана, но исследование по обоснованности применения этого критерия уже было проведено. А до этого исследования, если не считать случая по оптимизации коэффициентов в формуле Планка, я оптимизировал параметры систем только по отклику системы, а не по разнице между откликом и каким-то заданным  значением. Но теперь получается, что в предыдущей статье я уже исследовал на практике критерий (4-4), когда мы вели оптимизацию по разнице между расчетными координатами, т.е. полученными при вычислительном эксперименте YRas(U) и наблюдаемыми координатами YNab(U) в каждом U-ом эксперименте плана. Этот критерий я называю dY в противовес критерию Y, который обычно применяется при многофакторном планировании, и где, при проведении натурных экспериментов, Yu0(U)= YNab(U), а, при проведение вычислительных экспериментов, Yu0(U)= YRas(U).

 

 

Но сложность поставленной перед нами задачи не ограничивается только заменой критерия Y на dY. Мы имеем многокритериальную задачу, где аномальный остаток по смещению каждого параметра является самостоятельным критерием оптимизации, а в таком случае оптимальные параметры по одному критерию могут оказаться совсем не оптимальными по другому. Вспомните Леверье, который для теории каждой из четырех планет получал разные оптимальные значения масс других планет. И чаще всего нам и приходится решать многокритериальные задачи, где свести все частные критерии к одному никак не получается. Скажу больше, получение единых критериев оптимизации при решение многокритериальных задач вопрос даже более сложный, чем проблема вековых смещений параметров орбит, и я обычно терпеть не могу когда используют единые критерии оптимизации, полученные по принципу автомотовелофотопримусрадиоружье или с применением весовых коэффициентов, т.е. или когда один критерий умножают на другой, а потом делят на третий  и т.д. или выставляют экспертные оценки (весовые коэффициенты) различным по природе критериям. Правда подобными фокусами, где процветает махровый субъективизм, чаще всего пользуются исследователи социально-экономических или технико-экономических систем. И мне даже пришлось коренным образом модернизировать функционально-стоимостной анализ, чтобы можно было объективным путем получить единый критерий оптимизации для социально-экономических или технико-экономических систем. Но в нашем случае нам идеально подойдет формула (4-1), приведенная в работах [12] и [13] для расчета комплексного (единого) критерия оптимизации, т.к. мы имеем дело с единой природой частных критериев оптимизации.

 

Хотя эта формула и не лишена недостатков. Например, для таких показателей как себестоимость, энергоемкость, металлоемкость и т.д. возможен абсурдный результат при расчете единого критерия оптимизации. Например, если в каком то эксперименте получится себестоимость 3 рубля, а оптимальное значение принято 4 рубля, то значение единого критерия в этом эксперименте получится хуже, чем в эксперименте, где себестоимость была 4 рубля. Но для оптимизации скорости распространения гравитации по вековым смещениям параметров орбит планет эта формула вполне подходит и в обоих экспериментах, где вековое смещение в вычислительном эксперименте окажется как больше наблюдаемого, так и меньше, значение критерия оптимизация окажется хуже, чем в эксперименте, где расчетное значение будет равно наблюдаемому. По этому, я буду использовать формулу (4-2) при проведении экспериментов по поиску скорости гравитации, а упрощенный вариант этой формулы (4-3) я буду использовать для исследования применимости критерия dY. При этом использовать как в формуле (4-4) абсолютные отклонения наблюдаемых и расчетных значений нельзя, т.к. абсолютные значения вековых смещений у разных параметров не сопоставимы по величине и имеют даже разную размерность. А весовые коэффициенты, при этом, будут у меня отражать достоверность аномальных отклонений наблюдаемых смещений параметров орбит от расчетных и чем больше достоверность этих отклонений тем выше я буду давать вес по этому параметру для этой планеты. А, например, по перигелию Венеры я приму вес равный нулю, т.к. более-менее достоверно определить наблюдаемое смещение перигелия Венеры нам вообще не удалось (также как и по восходящему узлу Земли в интервале с 1800 по 2000 год).

 

Yu0(U) = SUMi ( kVesa(I) * ((YRas(I, U) - Yopt(I)) / Yopt(I)))^2                          (4-1)

где: Yu0(U) – разница между расчетным или наблюдаемым значением критерия в U-ом эксперименте для всех I–х параметров

Yopt(I) – оптимальное или рекомендуемое значение I-ого параметра

kVesa(I) – весовой коэффициент I-ого параметра (от 0 до 1).

 

Yu0(U) = SUMi,j ( kVesa(I, J) * Abs((YRas(I, J, U) - YNab(I, J)) / YNab(I, J)) / 100)     (4-2)

где: Yu0(U) – сумма разниц между расчетными YRas(I, J, U) и наблюдаемыми YNab(I, J) значениями вековых смещений в U-ом эксперименте для J –ой планеты и I -го параметра

kVesa(I, J) – весовой коэффициент I -ого параметра для J–ой планеты (от 0% до 100%).

 

Yu0(U) = Abs((YRas(U) - Yopt) / Yopt)                                         (4-3)

где: Yu0(U) – относительная разница между расчетным YRas(U) и оптимальным Yopt откликом системы, рассчитанным по одной из формул (3t…10t), которые имитируют поведение какой-то системы в U-ом эксперименте.

 

Yu0(U) = (SUMi Abs((YRas(U,i) - YNab(U,i))) / Ni                    (4-4)

где: YRas(U,i) =R(X1(U,i),Y1(U,i),Z1(U,i)) – расчетные координаты i-го наблюдения

        YNab(U,i) =R(X2(U,i),Y2(U,i),Z2(U,i)) – наблюдаемые координаты i-го наблюдения

        i – номер наблюдения в выборке из Ni наблюдений

 

Но, если вопросы многокритериальности хотя бы обсуждаются в теории оптимизации, то вопрос о том, какая аппроксимация по отклику системы или по отклонению отклика от заданного значения более эффективно и точно позволяет найти оптимальные параметры, в работах по многофакторному планированию вообще не освящен. А мой собственный опыт с коэффициентами в формуле Планка (формула 1t) не очень показательный, т.к. оптимизировал я там только 3 коэффициента, а 4-ый фактор (температура излучения) мною принудительно задавался для повышения качества информации, при проведении вычислительных экспериментов по почти D-оптимальному плану Бокса для четырех факторов. К тому же второй коэффициент (показатель степени при частоте излучения - v) мне надо было не столько оптимизировать, сколько подтвердить, что он равен 3, как это следовало из формулы Вина. Таким образом, я по большому счету оптимизировал только 2 параметра (фактора) и, следовательно, у меня могли быть только двойные смешанные взаимодействия, а это прекрасно воспроизводится полиномом (уравнением регрессии) 2-ой степени (3), хотя в формуле (1t) присутствует и не квадратичный член.

 

Да, теория многофакторного планирования позволяет значительно снизить необходимое количество экспериментов, которые надо выполнить, чтобы получить квадратичное уравнение регрессии аналогичное уравнению (3). При этом априори принимается, что в исследуемой системе нет взаимодействий между факторами, влияющими на критерий оптимизации, выше парных, т.е. нет, например, взаимодействия X1*X2*X3 или X1*X2*X3*X4. По этому, приступая к исследованию какой то системы, желательно на тестовых примерах убедится, что в системе нет взаимодействий выше парных или их влияние на критерий оптимизации не велико. Но как себя практически проявит критерий dY, если в исследуемой системе есть такие взаимодействия не известно. По этому, я решил проверить точность оптимизации параметров по критериям dY и dY^2  на простейших математических выражениях (3t…10t), которые будут имитировать поведение различных систем. А в качестве критерия оптимизации я буду использовать выражение (4-3). Только в нем мною используется значение разности между откликом системы и оптимальным значением по абсолютной величине, и это отличается от формулы (4-1), но, как пишут авторы [12], это делается только для того, чтобы разность была всегда положительной. А, если мы возведем в квадрат выражение по формуле (4-3), то получим формулу (4-1) по этому, будет возможность исследовать точность оптимизации как по критерию dY так и по критерию dY^2.

 

YRas(U) = X1*v^X2* exp(-X3*v/X4) (1t)

YRas(U) = u (2t)

YRas(U) = X1 * X2 * X3 / X4^2 (3t)

YRas(U) = X1^2 * X2^2 (4t)

YRas(U) = X1 * X2 * X3 * X4 (5t)

YRas(U) = X1 + X1^2 (6t)

YRas(U) = X1 + X2 + X3 + X4 (7t)

YRas(U)  = X1*X2 + X1*X3 + X1*X4 + X2*X3 + X2*X4 + X3*X4 (8t)

YRas(U)  = X1^2 + X2^2 + X3^2 + X4^2 (9t)

YRas(U) = X1 + X2 + X3 + X4 + X1^2 + X2^2 + X3^2 + X4^2 +

+ X1*X2 + X1*X3 + X1*X4 + X2*X3 + X2*X4 + X3*X4 (10t)

 

И первым делом я решил взять чуть ли не самый сложный случай с четверным взаимодействием, где вдобавок одно взаимодействие еще и не линейно, т.е. со всем Вам известным законом тяготения Ньютона (3t) и попробовать оптимизировать его параметры с использованием почти D-оптимального плана Бокса для четырех факторов (этот план отличается от приведенного в предыдущей статье ортогонального плана отсутствием нулевой точки, т.е. 25-го эксперимента).  В принципе, мы можем с законом тяготения провести и натурные эксперименты. Или, не с самим законом тяготения для масс, а с законом тяготения для зарядов (закон Кулона), где даже аналог гравитационной постоянной можем изменять, располагая различные диэлектрики между зарядами. Но речь сейчас идет не о том, можем ли мы воспроизвести эксперименты на реальном объекте или на его модели, а о том, можем ли мы, уже даже зная аналитическую формулу, отражающую отклик системы на наши воздействия на нее, чисто с математической точки зрения получить оптимальные значения параметров системы по примененному мною критерию dY, т.е. по разнице между откликом системы и известным оптимальным значением.

 

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

 

Для исследования возможностей методов многофакторного планирования при использование различных критериев оптимизации я сделал специальную форму 18 (см. рис. 8), где можно это сделать на нескольких тестовых примерах, анологичных 1t...10t, если в рамке //экспериментальные данные// отметить чекбокс //доп. тестов//. А на рис. 8 Вы видите результат аппроксимации стандартного теста (синие и красные кружки), который будет аналогичен тесту 10t для заданного Вами количества факторов при использование критерия Y и схему решения линейных уравнений при оптимизации параметров системы. При этом, на рисунке мы видим только часть информации, а если прокрутить рисунок, то Вы увидите и план эксперимента и получившиеся коэффициенты для уравнения (3). Эту форму программы Вы можете использовать и для обработки произвольных данных полученных вами при проведении каких то своих вычислительных или натурных экспериментов. Для этого Вам надо создать план, выбрав вид плана и количество факторов, и выполнить по этому плану эксперименты, а затем записать полученные результаты в файл. И перед тем как Вы нажмете кнопку //Выполнить эксперименты и рассчитать критерий// надо будет в рамке //экспериментальные данные// отметить переключатель //из файла//.

 

 

Рис.8 Исследование возможностей различных планов и по разным критериям оптимизации. Скриншот программы Solsys7.

 

 

Результаты оптимизации параметров в формуле тяготения Ньютона (3t), по примененным мною критериям оптимизации Y (нижняя часть рисунка) и dY (верхняя часть) представлены на рис. 9. Здесь маленькие кружки это критерий оптимизации, рассчитанный с использованием имитатора (3t) по формуле (4-3) в 24 экспериментах плана Бокса (номер соответствует абсциссе), а большие  кружки это этот же критерий оптимизации, но рассчитанный по уравнению регрессии (3), для тех же значений параметров, что и в соответствующем эксперименте.  При этом все факторы X1 … X4, при выполнение плана Бокса, на нулевом уровне были равны единице, а интервалы их варьирования в различных вариантах расчета были 0,2, 0,5 и 0,8 и значение Yopt бралось равным 1, т.е. я принимал, что оптимальные значения параметров X1 … X4 в формуле (3t) равны 1 и получалось, что Yopt = 1*1*1/1=1.

 

 

Рис.9 Точность аппроксимации закона притяжения Ньютона по разным критериям.

 

Как видим, по критерию оптимизации Y, т.е. по самому отклику системы, при небольших интервалах варьирования результат получается хороший, не смотря на то, что уравнение (3t) не должно аппроксимироваться уравнением (3), но вот при средних и больших интервалах варьирования результаты получаются не очень впечатляющие (чего и стоило ожидать), а по критерию dY результаты получаются еще хуже. По этому я решил продолжить исследование с использованием других математических выражений имитирующих поведение системы. При проведении вычислительных экспериментов по всем этим выражениям я принимал, что оптимальные значения всех параметров в этих формулах равны единице и находил сначала оптимальное значение отклика системы Yopt, а затем задавал в соответствие с планом различные значения параметров X1 … X4 и, вычислив значение YRas(U), находил значение критерия оптимизации в U-ом эксперименте. Значение критерия Y я брал равным YRas(U), критерия dY я вычислял по формуле (4-3), а критерий dY^2  определял, возведя критерий dY в квадрат. Затем, обработав данные по 24 значениям этих критериев Yu0(U), я получал коэффициенты k0 … k14 для аппроксимации (3) по которой вычислял значение Yu(U) для тех же значений параметров X1 … X4, что были в 24 экспериментах по плану эксперимента, и сравнивал их с 24 значениями Yu0(U). А по результатам сравнения я выставлял оценку почти D-оптимальному плану Бокса по качеству аппроксимации экспериментальных данных выражением (3) по различным критериям оптимизации.

 

Сравнение качества аппроксимации я проводил графически определяя попало ли значение Yu0(U), которое на приведенных выше рисунках отражено в 24 экспериментах маленькими кружками, внутрь большого кружка отражающего значение Yu(U). При этом, графический масштаб для вывода данных выбирался так, чтобы все данные от минимального до максимального значения по ординате укладывались в интервале от 1 до 7 сантиметров. А оценки качеству аппроксимации, т.е.  полученным по плану Бокса полиномам (3) при использовании критериев оптимизации Y, dY и dY^2 , я производил по пятибалльной шкале, определяя наилучшую оценку по наихудшим результатам, а затем, полученные результаты, оформил в виде таблицы. Критерии оценок были такими

5 баллов – все 24 маленьких кружка, т.е. значения Yu0(U), находятся внутри соответствующего большого кружка, центр которого соответствует ординате значения Yu(U). При этом диаметр большого кружка в два раза больше диаметра маленького.

4 балла - все 24 маленьких кружка или находятся внутри больших кружков или хотя бы касаются его с наружной стороны. При этом я также указываю в таблице, рядом с оценкой, в знаменателе количество экспериментов с наихудшими результатами. Так, если в 23 случаях маленькие кружки находятся внутри больших (это 5 баллов), а в одном случае маленький кружок пересекается с большим или касается его с наружной стороны (это 4 балла), то будет указана оценка 4/1.  

3 балла – если все маленькие кружки находятся хотя бы на расстояние 1-го диаметра большого кружка от его окружности. При этом, если 20 маленьких кружков находятся внутри больших, 1 кружок пересекается с большой окружностью, а три маленьких кружка находятся на расстояние 1-го большого диаметра от его окружности, то оценка будет 3/3. Таким образом, в знаменателе указывается только количество кружков не прошедших по более высоким оценкам.

2 балла – если все маленькие кружки находятся хотя бы на расстояние 2-х диаметров большого кружка от его окружности.

1 балл - если хотя бы один маленький кружок находится на расстояние более 2-х диаметров большого кружка от его окружности.

 

Таблица 12. Оценки качества аппроксимации полиномом 2-ой степени различных имитаторов  (имитирующих поведение системы) по различным критериям оптимизации при проведении вычислительных экспериментов по почти D-оптимальному плану Бокса при разных интервалах варьирования оптимизируемых факторов (параметров), которые на нулевых уровнях принимали значение 1. 

 

формула

критерий Y

критерий dY

критерий dY^2

 

+/-0,2

+/-0,5

+/-0,8

+/-0,2

+/-0,5

+/-0,8

+/-0,2

+/-0,5

+/-0,8

3t

4 / 2

2 / 1

1 / 2

3 / 3

2 / 2

1 / 2

1 / 2

1 / 3

1 / 6

4t

5

4 / 4

3 / 4

3 / 4

3 / 4

3 / 4

2 / 4

2 / 4

2 / 4

5t

5

3 / 10

1 / 1

2 / 1

1 / 2

1 / 6

1 / 2

1 / 6

1 / 6

6t

5

5

5

5

5

5

5

5

5

7t

5

5

5

3 / 16

3 / 16

3 / 16

5

5

5

8t

5

5

5

2 / 1

2 / 1

2 / 1

3 / 2

2 / 2

1 / 2

9t

5

5

5

3 / 11

3 / 5

3 / 6

4 / 8

4 / 9

3 / 6

10t

5

5

5

3 / 16

3 / 15

3 / 15

4 /2

3 / 2

2 / 2

 

Что можно сказать о полученных данных. То, что в тестах 6t…10t по критерию Y результаты будут отличными я и не сомневался, т.к. для аппроксимации именно таких взаимодействий полиномы 2-го порядка и приспособлены, но вот почему по критериям dY и dY^2 получились совсем плачевные результаты это для меня неожиданность. Ясно было и то, что в тестах 3t…5t результаты будут хуже чем в тестах 6t…10t, т.к. линейные взаимодействия выше парных в примененном нами плане Бокса смешаны с другими взаимодействиями и по этому их нельзя выделить (про парное квадратичное взаимодействие вообще затрудняюсь сказать что то определенное), но и в этом случае я никак не ожидал, что по критериям dY и dY^2 будут такие плохие результаты. Причем по результатам тестирования затруднительно даже сказать какой из этих критериев более точно позволяет найти оптимум, т.к. в одних тестах лучшие результаты показал критерий dY, а в других dY^2.  По этому остается только надеяться на то, что в исследуемой нами Солнечной системе тройных и выше взаимодействий не будет или их влияние будет не значительно, а то нам придется очень долго мучиться пока мы найдем область оптимума.

 

 

Однако, давайте наконец то посмотрим насколько удачно наш полином (3) аппроксимирует смоделированные нами процессы в Солнечной системе по критериям dY и dY^2. На рисунке 10 (в верхней части) показаны значения Yu0(U) полученные после обработки данных вычислительных экспериментов на классической математической модели Солнечной системы с учетом скорости распространения гравитации по критерию dY по формуле (4-2) (маленькие синие кружки) и по критерию dY^2 (маленькие красные кружки) и соответствующие им значения Yu(U) (большие кружки), полученные по уравнениям (3). А если на рисунке Вы видите маленький или большой кружок только одного цвета, значит их ординаты совпали. Причем в уравнение (4-2) значения рассчитаны по комплексному критерию с учетом вековых смещений всех I -х параметров, но от одного Меркурия, т.е. I =1 … 4, а J =1 и с весовыми коэффициентами kVesa(I, J)=100. Как видим по критерию dY аппроксимацию можно оценить на твердую тройку, а если бы не 23 и 24 эксперименты, где у нас был очень большой интервал варьирования скорости гравитации, то возможно бы получилась и четверка. Кстати и интервал варьирования скорости системы по оси Z (21 и 22 эксперименты) тоже великоват, т.к. в экспериментах 17 и 18, 19 и 20, 21 и 22, 23 и 24 разность значений критерия между этими парными экспериментами должна быть примерно равна. И по критерию dY^2 аппроксимацию тоже можно оценить на тройку, а если уменьшить интервалы варьирования скорости системы по оси Z и скорости гравитации то тоже может получиться хорошая оценка. Таким образом, быстрее всего в исследуемой нами Солнечной системе взаимодействия выше парных есть, но их влияние на критерий оптимизации все же мало, и по этому для оптимизации скорости гравитации и абсолютной скорости Солнечной системы вполне возможно использовать планы многофакторного планирования второго порядка, как с критерием dY, так и с критерием dY^2, но я буду применять критерий dY.

 

 

 

Рис.10 Оценка точности аппроксимации процессов в Солнечной системе и в тесте 2t.

 

Кстати, если кому-то показалось (глядя на выставленные мною оценки по тестам), что методы многофакторного планирования не очень надежны, то хочу Вас разочаровать, т.к. в основе этих методов лежит очень мощный математический аппарат, хотя и зародились эти методы для исследований в сельском хозяйстве, т.е. в области не очень блещущей математическими изысками. А в доказательство хочу Вам показать маленький фокус, на котором я частенько тестирую свои новые программы, где применяю многофакторное планирование. На нижней части рисунка 10 Вы видите значения Yu0(U) по критерию Y для теста 2t и эти же значения, полученные по уравнению регрессии (3) для Yu(U). А фокус заключается в том, что не зависимо от уровней варьирования факторов мы принимаем значение Yu0(U) в очередном эксперименте просто равным номеру этого эксперимента, т.е. получается, что уравнением (3) мы аппроксимировали почти полный хаос, но результат то получился вполне приличным. А, если мы найдем коэффициенты в формуле (3) по тестовому уравнению 10t по критерию оптимизации Y, то мы получим просто отличный результат (см. рис 8) и коэффициенты k0...k14 у нас получатся равны 1 с точностью до 7-го знака, т.е. претензий к самим методам многофакторного планирования никаких быть не может и если у нас и возникнут какие-то проблемы при поиске скорости гравитации, то их причина не методы многофакторрного планирования.

 

 

 

 

                                                                 Список литературы

 

1. – В.А.Брумберг //Релятивистская небесная механика// М.: Наука, 1972, 382 с.

2. – Кеплеровы элементы орбит с 1800 по 2050годы  http://ssd.jpl.nasa.gov/txt/p_elem_t1.txt

3.- Г.А.Чеботарев //Аналитические и численные методы небесной механики// М.: Наука, 1965, 367 с.

4.- С.Ю.Юдин //Кинематическая теория планет//  2012, 72 с http://modsys.narod.ru/Stat/Stat_Est/Kinematik/Kinematik1.zip

5. – С.Ю.Юдин //Модели и имитаторы// труды 5-ой международной научно-технической конференции //Компьютерное моделирование 2004// Часть 1. СПб.: Нестор, 2004, 356 с. http://modsys.narod.ru/Stat/Stat_Est/Konf_SPB/mod_imi.html

6. – S. Newcomb //Four inner planets and the fundamental constans of astronomy// - Washington, Goverment printing office, 1895, 202 p.

7. – М.Ф.Субботин //Введение в теоретическую астрономию// - М.: Наука, 1968, 800 с.

8. – Карим Хайдаров //Эфир великий часовщик// http://www.bourabai.kz/watchmaker.htm

9. – Н.Т.Роузвер //Перигелий Меркурия. От Леверье до Эйнштейна//: пер. с англ. – М.: Мир, 1985. – 246 с.

10. – О.В.Зайцев //Принцип Маха и орбитальная прецессия планет //                           http://www.physical-congress.spb.ru/russian/pmoop/pmoop.asp#_Toc479550256

11. - Справочное руководство по небесной механике и астродинамике, под редакцией Г.Н. Дубошина, издание второе, дополненное и переработанное - М.: Наука, 1976. – 864 с.

12. -  С.В.Мельников, В.Р.Алешкин, П.М.Рощин //Планирование эксперимента в исследованиях сельскохозяйственных процессов// Л.: Колос, 1980, 168 с.

13 – Адлер Ю.П., Маркова Е.В., Грановский Ю.Б. // Планирование экспериментов при поиске оптимальных условий// М.: Наука, 1967, 280 с.

14. –  В.Л.Гинзбург //Экспериментальная проверка общей теории относительности// Успехи физических наук, т. LIX, вып. I, Май, 1956, 47 с.

15. – С.Ю.Юдин //Моделирование систем и оптимизация их параметров// - Волгоград: Электронный вариант книги, 2003. - 208с.  http://modsys.narod.ru/Stat/kniga.html

 

 

Hosted by uCoz