Численное моделирование и выбор шага.

[akkoch]

Добрый день, друзья!

Помогите разобраться с проблемой.
Рисую небольшую программку с ребенком. Модель условно Солнечной системы.
Уравнения движения решаем методом Эйлера (да, метод не точный, но для начала этого хватало)
Решили заменить метод Эйлера на метод Рунге-Кутты 4 порядка.
В нашей модели есть "центральное тело" - некая звезда массы M планета массы m и спутник планеты.
Когда решали методом Эйлера, высокая скорость планеты (кратно большая ее димаметру) не как не мешала:
1) запоминали положение планеты S0
2) двигали планету относительно центрального тела
3) двигали спутник относительно S0 планеты.
Все расчеты ведем в векторном виде.
Проблема возникала при попытке перехода на РК4. Решили внедрить его для движения спутника.
Мы как и раньше:
1) запоминали положение планеты S0
2) двигали планету относительно центрального тела
3) при вычислении k1...k4 мы считаем ускорение в этих промежуточных точках, и из-за высокой скорости эти промежуточные точки оказываются на значительном удалении от S0 и движение спутника превращается в кошмар :)
Я думаю что проблема в этом, при ближе скорость движения планеты к нулю, тем ближе решение к решению методом Эйлера, либо надо шаг приращения уменьшать на несколько порядков.
Мы хотели при использовании более точного метода увеличить шаг интегрирования, но получили облатный эффект :)

Хочу вашего совета. Как в такой ситуации поступить?
Если при вычислении k1..k4 для спутника, считать гравитацию не от S0, а от каких то S0+k1, S0+k2... для планеты, те считать эти коэффизиенты для спутника и планеты условно параллельно.
Или в подобных моделях используется совершенно другой подход?

Спасибо!

[Pphantom]

akkoch в сообщении #1491052 писал(а):

В нашей модели есть "центральное тело" - некая звезда массы M планета массы m и спутник планеты.

А какими были реальные использованные величины? И массы, и начальные координаты и скорости?

А то по описанию возникают нехорошие подозрения, что вы попросту решаете не ту задачу (и "кошмар" во втором случае больше похож на правильное решение, чем "ожидаемое" поведение - в первом).

[Pphantom]

Тема перемещена из форума «Помогите решить / разобраться (М)» в форум «Астрономия»
Причина переноса: и вообще поехали в профильный раздел.

[B@R5uk]

akkoch, а по-подробнее формулы, по которым вы считаете, можете изложить? А то не понятно, что с чем сравнивать и что проверять.

На всякий случай прокомментирую сам метод Рунге-Кутта. Он решает задачу

$$\frac{d\mathbf{R}}{dt}=\mathbf{F}\left(\mathbf{R},t\right)$$

Здесь $\mathbf{R}$ — столбец, который содержит и координаты, и скорости всех движущихся тел. Функция $\mathbf{F}$ — тоже столбец, который содержит, соответственно, скорости и ускорения. Так что коэффициенты $k_1,\,...,\,k_4$ из метода Рунге-Кутта так же являются столбцами и содержат скорости и ускорения.

[Утундрий]

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

[Emergency]

akkoch в сообщении #1491052 писал(а):

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

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

[B@R5uk]

Emergency в сообщении #1491089 писал(а):

Тут сэкономить не получится.

Вот те раз! А зачем тогда эти численные методы напридумывали, если не из-за экономии? Метод Эйлера вычисляет функцию один раз за шаг, но является методом 1-го порядка. Чтобы увеличить точность в 100 раз необходимо производить в 100 раз больше обращений к функции. Метод Рунге-Кутта 4-го порядка является методом 4-го порядка (извините за тавтологию). За шаг он вычисляет функцию 4 раза, за то, чтобы увеличить точность в 100 раз необходимо сократить шаг всего в $\sqrt{10}\approx 3,2$ раза. Разве не экономия? Я, конечно, понимаю, что всего лишь порядки, и при переходе от одного метода к другому нужно учесть абсолютный коэффициент потери/набора точности, но тем не менее.

-- 07.11.2020, 21:30 --

Тут другое дело. Если метод Эйлера не сильно чувствителен к ошибкам, то метод Рунге-Кутта ошибок не прощает. Другими словами, если в методе Эйлера сначала прибавлять к скорости, а потом координате, а не наоборот, как положено, то точность не пострадает, в то время как при наличии подобной ошибки в методе Рунге-Кутта, порядок метода упадёт с 4-го сразу до 1-го. Все эти лишние расчёты уйдут псу под хвост.

[Emergency]

B@R5uk в сообщении #1491099 писал(а):

Вот те раз! А зачем тогда эти численные методы напридумывали, если не из-за экономии?

Наверное я не совсем правильно понял задачу. Я понял ее так, что ТС получает результат в виде отображения траектории на мониторе. Если монитор высокого разрешения, размер траектории по вертикали можно принять за 1000 писелей. Тогда ее длина будет примерно 3000 пикселей, координаты каждого из которых требуется вычислить. На это потребуются вычислительные ресурсы, про которые я и говорил.

[DimaM]

akkoch
Для контроля шага хорошо бы следить за интегралами движения: суммарной энергией и суммарным моментом импульса. Они должны сохраняться с точностью, скажем, $10^{-3}$. Если расхождение больше, шаг нужно уменьшать.
RK4 для вашей задачи видится мне перебором. Формул много, легко где-нибудь ошибиться.
Я бы посоветовал Эйлера-Кромера (координаты изменяем, используя начальные скорости; затем изменяем скорости, используя ускорения, рассчитанные по измененным координатам). Метод хорош тем, что точно сохраняет момент импульса.
Если хочется метода повышенной точности, можно взять Верле в скоростной форме. Он имеет второй порядок, но требует вычисления силы только один раз за шаг.

[akkoch]

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

Pphantom в сообщении #1491061 писал(а):

А какими были реальные использованные величины? И массы, и начальные координаты и скорости?

А то по описанию возникают нехорошие подозрения, что вы попросту решаете не ту задачу (и "кошмар" во втором случае больше похож на правильное решение, чем "ожидаемое" поведение - в первом).

Центральное тело имеет массу Солнца. "планета" - тело 30км в радиусе, плотностью $600 kg/m^{3}$. Что-то вроде кометного вещества.
Уточню на счет ожидаемого результата: пробовали просто "подбрасывать" некую частицу с поверхности (считая методом Эйлера) семитрично с разных полусфер направляя "планету" по прамой на центральное тело ("планета" сферическое тело). И считая указанным методом мы получаем семметричные траектории подброшенных частиц. Что ожидаемо, и косвенно указывает на правильность расчетов. Более того такой же симметричный результат, получам и методом более высокого порядка, но только если начальную скорость планеты положить за ноль.

DimaM в сообщении #1491202 писал(а):

Я бы посоветовал Эйлера-Кромера (координаты изменяем, используя начальные скорости; затем изменяем скорости, используя ускорения, рассчитанные по измененным координатам). Метод хорош тем, что точно сохраняет момент импульса.
Если хочется метода повышенной точности, можно взять Верле в скоростной форме. Он имеет второй порядок, но требует вычисления силы только один раз за шаг.

Спасибо за ваши рекомндации! Обязательно реализуем.

DimaM в сообщении #1491202 писал(а):

Для контроля шага хорошо бы следить за интегралами движения: суммарной энергией и суммарным моментом импульса. Они должны сохраняться с точностью, скажем, $10^{-3}$.

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

Emergency в сообщении #1491157 писал(а):

Я понял ее так, что ТС получает результат в виде отображения траектории на мониторе. Если монитор высокого разрешения, размер траектории по вертикали можно принять за 1000 писелей.

Вообще с графикой не заморачиваемся пока. Результат получаем в текстовый файл, который потом визуализируем (сначала был Excel, теперь Grapher)

Утундрий в сообщении #1491076 писал(а):

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

Да всё просто - увлечены астрономией, а что-то сделать и(или) разобраться это на пользу подрастающему поколению.
я тут вообще инкогнито :) Конечно, если идеи и советы полученные от вас помогут нам продвитуться, я скажу что это благодаря вам! Давным давно еще в студенчестве я и кодил и матан изучал, потом конечно ничем этим почти не пользовался. И вот до определенного уровня моих знаний пока нам хватало.

B@R5uk в сообщении #1491064 писал(а):

а по-подробнее формулы, по которым вы считаете, можете изложить?

Это конечно можно. Но здесь, если не ошибаюсь, в Латехе они должны быть. Тут нужно отдельно поработать, поэтому я могу пошагово описать алгоритм словами, если так можно :)

B@R5uk в сообщении #1491064 писал(а):

Так что коэффициенты $k_1,\,...,\,k_4$ из метода Рунге-Кутта так же являются столбцами и содержат скорости и ускорения.

Если я правильно понял и вас и метод, то коэффициенты $k_1,\,...,\,k_4$ мы потом интегрируем по времени и получем столбец координат и скоростей.


Еще раз, спасибо за отзывы!

[Pphantom]

akkoch в сообщении #1491369 писал(а):

Центральное тело имеет массу Солнца. "планета" - тело 30км в радиусе, плотностью $600 kg/m^{3}$. Что-то вроде кометного вещества.

А спутник? Начальное положение относительно планеты, начальная скорость?

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

В общем, опишите постановку задачи полностью - в т.ч. количественно.

[B@R5uk]

akkoch в сообщении #1491369 писал(а):

Если я правильно понял и вас и метод, то коэффициенты $k_1,\,...,\,k_4$ мы потом интегрируем по времени и получем столбец координат и скоростей.

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

akkoch в сообщении #1491369 писал(а):

Но здесь, если не ошибаюсь, в Латехе они должны быть.

Проще код программы вставить.

[B@R5uk]

Как-то так можно это сделать:

Используется синтаксис Java

import java.util.*;

class VectorClass {
        
        private double [] data;
        
        VectorClass () {
                data = new double [3];
        }
        
        VectorClass (double [] arg) {
                data = Arrays .copyOf (arg, 3);
        }
        
        VectorClass (VectorClass arg) {
                data = Arrays .copyOf (arg .data, 3);
        }
        
        public double getModuleSquare () {
                double result = 0.0;
                for (int k = 0; 3 > k; ++k) {
                        result += data [k] * data [k];
                }
                return result;
        }
        
        public double getModule () {
                return Math .sqrt (getModuleSquare ());
        }
        
        public double getModuleCube () {
                double value = getModuleSquare ();
                return value * Math .sqrt (value);
        }
        
        public void add (VectorClass arg) {
                for (int k = 0; 3 > k; ++k) {
                        data [k] += arg .data [k];
                }
        }
        
        public void addWithConstant (double arg1, VectorClass arg2) {
                for (int k = 0; 3 > k; ++k) {
                        data [k] += arg1 * arg2 .data [k];
                }
        }
        
        public void multiplyBy (double arg) {
                for (int k = 0; 3 > k; ++k) {
                        data [k] *= arg;
                }
        }
        
        public VectorClass getVectorTo (VectorClass arg) {
                VectorClass result = new VectorClass (arg);
                for (int k = 0; 3 > k; ++k) {
                        result .data [k] -= data [k];
                }
                return result;
        }
        
        public double [] getCoordinates () {
                return Arrays .copyOf (data, 3);
        }
}

class BodyClass {
        
        public static final double graviConstant = 6.67430E-11;
        public static final double solarMass = 1.98847E30;
        public static final double earthMass = 5.9722E24;
        public static final double lunaMass = 7.342E22;
        public static final double au = 1.495978707E11;
        public static final double lunaOrbit = 3.844E8;
        
        private double mass;
        private VectorClass coordinates;
        private VectorClass velocity;
        
        BodyClass (BodyClass arg) {
                mass = arg .mass;
                coordinates = new VectorClass (arg .coordinates);
                velocity = new VectorClass (arg .velocity);
        }
        
        BodyClass (double arg) {
                mass = arg;
                coordinates = new VectorClass ();
                velocity = new VectorClass ();
        }
        
        BodyClass (double arg1, VectorClass arg2) {
                mass = arg1;
                coordinates = new VectorClass (arg2);
                velocity = new VectorClass ();
        }
        
        BodyClass (double arg1, VectorClass arg2, VectorClass arg3) {
                mass = arg1;
                coordinates = new VectorClass (arg2);
                velocity = new VectorClass (arg3);
        }
        
        public VectorClass getAccelFrom (BodyClass arg) {
                VectorClass result = coordinates .getVectorTo (arg .coordinates);
                result .multiplyBy (graviConstant * arg .mass / result .getModuleCube ());
                return result;
        }
        
        public VectorClass getAccelFrom (BodyClass [] arg) {
                VectorClass result = new VectorClass ();
                for (int k = 0; arg .length > k; ++k) {
                        if (this != arg [k]) {
                                result .add (getAccelFrom (arg [k]));
                        }
                }
                return result;
        }
        
        public VectorClass getVelocity () {
                return new VectorClass (velocity);
        }
        
        public void addWithConstant (double arg1, BodyClass arg2) {
                coordinates .addWithConstant (arg1, arg2 .coordinates);
                velocity .addWithConstant (arg1, arg2 .velocity);
        }
        
        public double [] getCoordinates () {
                return coordinates .getCoordinates ();
        }
}

class SystemClass {
        
        private BodyClass [] bodies;
        
        public SystemClass (SystemClass arg) {
                bodies = new BodyClass [arg .bodies .length];
                for (int k = 0; arg .bodies .length > k; ++k) {
                        bodies [k] = new BodyClass (arg .bodies [k]);
                }
        }
        
        public SystemClass (BodyClass [] arg) {
                bodies = Arrays .copyOf(arg, arg .length);
        }
        
        public SystemClass getFunction () {
                BodyClass [] newBodies = new BodyClass [bodies .length];
                for (int k = 0; bodies .length > k; ++k) {
                        newBodies [k] = new BodyClass (0.0, bodies [k] .getVelocity (), bodies [k] .getAccelFrom (bodies));
                }
                return new SystemClass (newBodies);
        }
        
        public void addWithConstant (double arg1, SystemClass arg2) {
                for (int k = 0; bodies .length > k; ++k) {
                        bodies [k] .addWithConstant (arg1, arg2 .bodies [k]);
                }
        }
        
        public double [] getCoordinates () {
                int m = 1;
                double [] result, temp;
                result = new double [1 + 3 * bodies .length];
                for (int k = 0; bodies .length > k; ++k) {
                        temp = bodies [k] .getCoordinates ();
                        for (int l = 0; 3 > l; ++l) {
                                result [m] = temp [l];
                                ++m;
                        }
                }
                return result;
        }
}

class RungeKuttaMethod {
        
        public static void update (SystemClass arg, double timeStep) {
                SystemClass K1, J2, K2, J3, K3, J4, K4;
                K1 = arg .getFunction ();
                J2 = new SystemClass (arg);
                J2 .addWithConstant (timeStep / 2.0, K1);
                K2 = J2 .getFunction ();
                J3 = new SystemClass (arg);
                J3 .addWithConstant (timeStep / 2.0, K2);
                K3 = J3 .getFunction ();
                J4 = new SystemClass (arg);
                J4 .addWithConstant (timeStep, K3);
                K4 = J4 .getFunction ();
                arg .addWithConstant (timeStep / 6.0, K1);
                arg .addWithConstant (timeStep / 3.0, K2);
                arg .addWithConstant (timeStep / 3.0, K3);
                arg .addWithConstant (timeStep / 6.0, K4);
        }
        
}

public class Simple_Planets_Modeling {
        
        private static int num = 2 * 370;
        private static final double [] [] data = new double [num] [];
        
        public static void main (String [] args) {
                double time, timeStep;
                SystemClass system = new SystemClass (new BodyClass [] {
                                new BodyClass (BodyClass .solarMass),
                                new BodyClass (
                                                BodyClass .earthMass,
                                                new VectorClass (new double [] {BodyClass .au, 0.0, 0.0}),
                                                new VectorClass (new double [] {0.0, 30000.0, 0.0})),
                                new BodyClass (
                                                BodyClass .lunaMass,
                                                new VectorClass (new double [] {BodyClass .au + BodyClass .lunaOrbit, 0.0, 0.0}),
                                                new VectorClass (new double [] {0.0, 31020.0, 0.0}))
                });
                time = 0.0;
                timeStep = 12.0 * 3600.0;
                for (int k = 0; num > k; ++k) {
                        time += timeStep;
                        RungeKuttaMethod .update (system, timeStep);
                        data [k] = system .getCoordinates ();
                        for (int l = 1; data [k] .length > l; ++l) {
                                data [k] [l] /= BodyClass .au;
                        }
                        data [k] [0] = time / (24.0 * 3600.0);
                        System .out .print (String .format ("%6.2f", data [k] [0]));
                        for (int l = 1; data [k] .length > l; ++l) {
                                System .out .print (String .format ("%11.6f", data [k] [l]));
                        }
                        System .out .println ();
                }
        }
}
 

Не рационально, конечно, куча классов постоянно по ходу вычислений создаётся и уничтожается, но за то структурировано и читабельно. И вроде работает. Причём метод Рунге-Кутта всего несколько строчек, всё остальное — операции над объектами. При аккуратном подходе любой метод будет работать, я только одну ошибку допустил во время написания: забыл самодействие объектов выключить. А так можно любое количество тел моделировать, коду всё равно сколько их при инициализации ему подсунешь.

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

[akkoch]

B@R5uk в сообщении #1491378 писал(а):

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

Ну значит я неверно выразился, конечно, для каждого шага мы берем в итоге $y_n_+_1 = y_n + \frac{h}{6} (k_1+2k_2+2k_3+k_4)$. И на следующем шаге $k_1$..$k_4$ вычисляем заново.

[Pphantom]

Тема перемещена из форума «Астрономия» в форум «Карантин»
по следующим причинам:

- поправьте все-таки формулы (краткие инструкции: «Краткий FAQ по тегу [math]» и видеоролик Как записывать формулы), и не только в последнем сообщении.

Исправьте все Ваши ошибки и сообщите об этом в теме Сообщение в карантине исправлено.
Настоятельно рекомендуется ознакомиться с темами Что такое карантин и что нужно делать, чтобы там оказаться и Правила научного форума.