Ограниченная симуляция 3 тел работает некорректно

В настоящее время я пытаюсь выполнить задание, где мне нужно написать симуляцию для ограниченной гравитационной задачи с тремя телами, с двумя фиксированными массами и одной тестовой массой. Информация, которую я получил по этой проблеме: Проверьте эту ссылку и вот моя программа на данный момент:

#include<stdlib.h>
#include<stdio.h>
#include <math.h>

int main (int argc, char* argv[])
{
    double dt=0.005, x[20000],y[20000],xv,yv,ax,ay,mneg,mpos,time,radius=0.01;
    int n,validation=0;
    FILE* output=fopen("proj1.out", "w");

    printf("\n");


    if((argv[1]==NULL) || (argv[2]==NULL) || (argv[3]==NULL) || (argv[4]==NULL) ||     (argv[5]==NULL) || (argv[6]==NULL)) 
    {
        printf("************************ ERROR! ***********************\n");
        printf("**     Not enough comand line arguments input.       **\n");
        printf("**     Please run again with the correct amount (6). **\n");
        printf("*******************************************************\n");
        validation=1;
        goto VALIDATIONFAIL;
    }
if((sscanf(argv[1], "%lf", &mneg)==NULL) || (sscanf(argv[2], "%lf", &mpos)==NULL) || (sscanf(argv[3], "%lf", &x[0])==NULL) ||
(sscanf(argv[4], "%lf", &y[0])==NULL) || (sscanf(argv[5], "%lf", &xv)==NULL) || (sscanf(argv[6], "%lf", &yv)==NULL) )
{
    printf("************************* ERROR! ************************\n");
    printf("** Input values must be numbers. Please run again with **\n");                 
    printf("** with numerical inputs (6).                          **\n");
    printf("*********************************************************\n");
    validation=1;
  goto VALIDATIONFAIL;
}

sscanf(argv[1], "%lf", &mneg);
sscanf(argv[2], "%lf", &mpos);
sscanf(argv[3], "%lf", &x[0]);
sscanf(argv[4], "%lf", &y[0]);
sscanf(argv[5], "%lf", &xv);
sscanf(argv[6], "%lf", &yv);


    x[1]=x[0]+(xv*dt);
    y[1]=y[0]+(yv*dt);

for(n=1;n<10000;n++)
    {
        if(x[n-1]>=(1-radius) && x[n-1]<=(1+radius) && y[n-1]>=(0-radius) && y[n-1]<=(0+radius))
        {
            printf("Test mass has collided with M+ at (1,0), Exiting...\n");
            goto EXIT;
        }
        else if(x[n-1]>=(-1-radius) && x[n-1]<=(-1+radius) && y[n-1]>=(0-radius) && y[n-1]<=(0+radius))
        {
            printf("Test mass has collided with M- at (-1,0), Exiting...\n");
            goto EXIT;
        }
        else
        {
            double dxn = x[n] + 1;
            double dxp = x[n] - 1;
            double mnegdist = pow((dxn*dxn + (y[n]*y[n])), -1.5);
            double mposdist = pow((dxp*dxp + (y[n]*y[n])), -1.5);

            ax =  -(mpos*dxp*mposdist+mneg*dxn*mnegdist);
            ay =  -(mpos*y[n]*mposdist+mneg*y[n]*mnegdist); 


            x[n+1]=((2*x[n])-x[n-1] +(dt*dt*ax));
            y[n+1]=((2*y[n])-y[n-1]+(dt*dt*ay));


            fprintf(output, "%lf %lf\n",x[n-1], y[n-1]); 
        }
    }
VALIDATIONFAIL: 
printf("\n");
return(EXIT_FAILURE);   

EXIT:
return(EXIT_SUCCESS);
}

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

Основная проблема заключается в том, что когда тестовая масса достигает точки на своей траектории, когда она должна оторваться и начать вращаться вокруг другой массы, вместо этого она просто улетает по прямой линии в бесконечность! Сначала я подумал, что массы сталкиваются, поэтому я включил проверку радиуса, но в некоторых случаях это работает, в некоторых нет, а в некоторых случаях массы сталкиваются раньше, прежде чем траектория все равно пойдет не так. так что дело явно не в этом. Я не уверен, что объяснил это слишком хорошо, поэтому вот картинка, чтобы показать вам, что я имею в виду. (моделирование справа взято из здесь)

строка

Однако это не всегда так, иногда вместо того, чтобы идти по прямой, траектория просто сходит с ума, когда она должна перейти к другой массе, вот так: сумасшедший

Я действительно абсолютно не знаю, что происходит, я потратил дни, пытаясь понять это, но просто не могу никуда добраться, поэтому любая помощь в определении моей проблемы будет очень признательна.


c simulation physics gravity orbit
person user1831711    schedule 19.12.2012    source источник
comment
Не имеет отношения к вашей проблеме, но для начала вы не должны сравнивать argv[x] с NULL, вместо этого проверьте argc. Кроме того, sscanf возвращает целое число, а не указатель, поэтому технически вы должны сравнивать с 0, а не с NULL. Вам также не нужно дважды вызывать sscanf для преобразования числа. И вместо использования sscanf для преобразования аргументов в значения с плавающей запятой проверьте strtod< /а> функция.   -  person Some programmer dude    schedule 19.12.2012
comment
Вероятно, есть куча вещей, которые вам нужно будет сделать, чтобы правильно кодировать, например, исправить отступ (это читается очень запутанно) и просто, пожалуйста, не используйте gotos. Возможно, ранние return могут быть плохими или запутанными, но goto уже давно ушли по пути Додо.   -  person    schedule 19.12.2012
comment
Тот факт, что ваша орбита совершенно сходит с ума, как только она попадает в один из центров тяжести, все равно указывает на деление на ноль. Я еще не проверил ваш радиус внимательно, но я думаю, что он все еще недостаточен. Еще один способ предотвратить деление на ноль — прибавить к знаменателю дроби очень минимальное значение, чтобы оно никогда не достигло нуля (здесь, конечно, это степень: x^-1,5, что равно (1/x ) ^ 1,5, так что все еще деление).   -  person    schedule 19.12.2012
comment
Вы пробовали с меньшим шагом времени? Когда ваша тестовая масса подходит очень близко к одной из звезд, потенциальная энергия преобразуется в кинетическую, и ее скорость становится очень высокой. При больших временных шагах это приводит к огромным ошибкам траектории, даже к явлениям выброса, которые вы наблюдаете. Это очень похоже на то, что мы наблюдаем в физике кластеров при моделировании молекулярной динамики со слишком большими временными шагами.   -  person Hristo Iliev    schedule 19.12.2012
comment
Христо, увеличение временного шага, кажется, остановило проблему выброса, однако траектория, которую я получаю, по-прежнему не совпадает с примером, который я опубликовал, хотя, по крайней мере, тестовая масса теперь переходит к массе 2.   -  person user1831711    schedule 19.12.2012
comment
и Эверт, добавление небольшого значения, как вы сказали, похоже, ничего не делает, к сожалению...   -  person user1831711    schedule 19.12.2012
comment
Задача трех тел хаотична, даже если два тела неподвижны; небольшое отклонение может быстро перерасти в большое. Теперь, когда @HristoIliev решил вашу проблему с выбросом, что заставляет вас думать, что ваши результаты неверны, а пример правильный?   -  person Beta    schedule 19.12.2012
comment
@Beta, на самом деле это хороший момент ... я думаю, я просто предположил, что, должно быть, я делаю это неправильно ... но теперь, глядя на это, мои результаты совпадают во всех случаях, кроме нескольких, с этим примером, так что я полагаю, что это может быть просто так они использовали немного другой метод или допустили небольшую ошибку. это неправильно ха-ха!   -  person user1831711    schedule 19.12.2012

Ответы (2)


arrow_upward
4
arrow_downward

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

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

Интегратор Верле, такой как реализованный вами, и схема Верле скорости имеют глобальную ошибку, которая составляет O(Δt2). Это означает, что алгоритм только квадратично чувствителен к временному шагу, и для значительного повышения точности необходимо соответственно уменьшить временной шаг во столько раз, сколько квадратный корень из желаемого коэффициента повышения точности. Нажмите кнопку в апплете Flash, с которым вы сравниваете свои траектории, и вы увидите, что он использует совершенно другой интегратор — алгоритм Эйлера-Кромера (также известный как полунеявный метод Эйлера). Она имеет разную точность при одном и том же временном шаге (на самом деле она хуже, чем у схемы Верле при одном и том же временном шаге), и, следовательно, вы не можете и не должны напрямую сравнивать обе траектории, а только их статистические свойства (например, средняя энергия, средняя скорость и т. д.)

Моя точка зрения заключалась в том, что вы должны уменьшить временной шаг, потому что он слишком велик для обработки случаев, когда тестовое тело подходит слишком близко к одному из центров гравитации. Здесь скрыта еще одна проблема, и это конечная числовая точность. Соблюдайте этот термин:

double dxp = x[n] - 1;
double mposdist = pow((dxp*dxp + (y[n]*y[n])), -1.5);

Всякий раз, когда вы вычитаете два близких числа с плавающей запятой (x[n] и 1.0), происходит очень неприятное событие, известное как потеря точности, поскольку большинство старших значащих битов в их мантиссе аннулируют друг друга, и в конце, после шага нормализации, вы получаете число с гораздо меньшими значащими битами, чем исходные два числа. Эта потеря точности становится еще больше, поскольку результат затем возводится в квадрат и используется в качестве знаменателя. Обратите внимание, что это в основном происходит вблизи оси симметрии системы, где y[n] приближается к 0. В противном случае y[n] может быть достаточно большим, так что dxp*dxp будет лишь небольшой поправкой к значению y[n]*y[n]. Конечным результатом является то, что сила будет совершенно неправильной вблизи каждой фиксированной массы и обычно будет больше по величине, чем фактическая сила. Это предотвращается в вашем случае, поскольку вы проверяете, находится ли точка за пределами предписанного radius.

Большие силы приводят к большим смещениям при заданном временном шаге. Это также приводит к искусственному увеличению полной энергии системы, т. е. пробная масса будет двигаться быстрее, чем при более тонком моделировании. Также может случиться так, что пробное тело окажется так близко к центру гравитации, что огромная сила, умноженная на квадрат временного шага, может привести к такому огромному смещению, что ваша пробная масса окажется намного дальше, но на этот раз с увеличение общей энергии приведет к высокой кинетической энергии, и тело будет практически выброшено из объема моделирования. Это также может произойти, даже если вы вычисляете силу с бесконечной точностью — просто смещение между двумя временными шагами может быть настолько большим (из-за большого временного шага), что система совершит нереалистичный скачок в фазовом пространстве к совершенно другой энергетической изоповерхности. А с гравитацией (как и с электростатикой) так легко дойти до такого случая, так как сила возрастает как 1/r^2 и вблизи радиуса на много порядков сильнее, чем в исходном состоянии.

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

Теперь, когда схемы Верле являются симплектическими, лучший способ контролировать правильность моделирования — это наблюдать за полной энергией системы. Обратите внимание, что интегратор скорости Верле немного лучше, поскольку он численно более стабилен (но все же имеет ту же зависимость точности от квадрата временного шага). Со стандартной схемой Верле вы можете получить приблизительную скорость v[i], взяв (x[i+1] - x[i-1])/(2*dt). При скорости Верле скорость явно включена в уравнения.

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

Уменьшение временного шага соответственно увеличивает время выполнения симуляции. Можно было также заметить, что в дальнем поле силы малы, точка движется медленно, а длинный временной шаг — это нормально. Более короткий временной шаг не улучшит решение, а только увеличит время выполнения. Вот почему люди изобрели алгоритмы с несколькими временными шагами, а также адаптивные алгоритмы временных шагов, которые автоматически уточняют решения в ближнем поле. Здесь также может быть применен другой метод вычисления сил путем преобразования уравнений таким образом, чтобы они не включали вычитание близких значений переменных.

(ну, это вышло намного больше, чем даже несколько комментариев)

person Hristo Iliev    schedule 19.12.2012
comment
Большое спасибо, вы даже не представляете, насколько вы были полезны. Вы не только решили мою проблему, но теперь я действительно понимаю ее и физику проблемы в целом. Ваша помощь была очень признательна, спасибо еще раз :) - person user1831711; 20.12.2012

arrow_upward
0
arrow_downward

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

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

F = ((-G * M1 * M2) / (R * R)) * r_unit_vector;

Где:

G ~= 6,67e-11, M1 — масса первого объекта, M2 — масса второго объекта, R — расстояние между двумя объектами: sqrt(pow(X2 - X1, 2) + pow(Y2 - Y1, 2))

X1 — координата X объекта 1, X2 — координата X объекта 2, Y1 — координата Y объекта 1, Y2 — координата Y объекта 2.

r_unit_vector — единичный вектор, указывающий от объекта 2 к объекту 1.

struct r_unit_vector_struct{
    double x, y;
}r_unit_vector;

r_unit_vector имеет компонент x, который является координатой x объекта 2 - координатой x объекта 1, r_unit_vector имеет компонент y, который является координатой y объекта 2 - координатой y объекта 1.

Чтобы сделать r_unit_vector единичным вектором, вы должны разделить (оба компонента x a и y отдельно) на его длину, которая определяется как sqrt(pow(r_unit_vector.x, 2) + pow(r_unit_vector.y - Y1, 2))

И вы все должны быть готовы к работе! Надеюсь, это имеет смысл. Если нет, я напишу вам класс, чтобы сделать это или объяснить это подробнее, если смогу!

person FreelanceConsultant    schedule 29.12.2012