расчет движения электронов в атоме лития методом МД(2мерная)

[Pphantom]

MathKvant в сообщении #974526 писал(а):

Спасибо, я вроде подправил самые тупые моменты,перевел все в СГСЭ,теперь коэф k=1,также вынес степень $10^{-10}$, теперь все компиллится, но значения все NaN.

Как-то странно Вы это сделали. У Вас вычисляется коэффициент const, но в дальнейшем он нигде не используется, вместо него во все выражения входит q*q/me. В чем смысл?

К тому же во все выражения заряды входят в квадрате. Частично можно было бы сказать, что это скомпенсировано использованием ангстремов в качестве единицы длины, но это $10^{-10}$ м, а не $10^{-10}$ см, у Вас остаются неучтенными два порядка.

А теперь внимательно посмотрите на начальные условия. В момент старта электроны №1 и №2 находятся в одной точке. Следовательно, при вычислении их ускорений Вы получите нуль в знаменателе. При этом для вычисления новых координат используются переменные вроде x_previous2, которые исходно все равны нулю.

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

[MathKvant]

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

Как-то странно Вы это сделали. У Вас вычисляется коэффициент const, но в дальнейшем он нигде не используется, вместо него во все выражения входит q*q/me. В чем смысл?

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

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

    program molecul
    implicit none
    real, parameter :: r1=0.52917720859! радиус первой орбиты в см*10^(-10)
    real, parameter :: r2=2.11670883436! радиус второй орбиты в см*10^(-10)
    real, parameter :: me=9.10938291e-18! масса электрона в гр*10^(-10)
    real, parameter :: q=4.8032042713! заряд элетрона в сгсэ
    real, parameter :: tau=1.52e-12! шаг по времени в секундах
    real v_next_x1,v_previous_x1,x_next1,x_previous1,v_next_y1,v_previous_y1,y_next1,y_previous1
    real v_next_x2,v_previous_x2,x_next2,x_previous2,v_next_y2,v_previous_y2,y_next2,y_previous2
    real v_next_x3,v_previous_x3,x_next3,x_previous3,v_next_y3,v_previous_y3,y_next3,y_previous3
    real x11,x22,x33,y11,y22,y33
    real, parameter ::  const=tau*q*q/me
    integer i
    write(*,*) const
    x11=r1! начальные условия
    x22=0
    x33=r2
    y11=0
    y22=r1
    y33=0
        v_previous_x1=0
        x_previous1=0
        v_previous_y1=0
        y_previous1=0


        v_previous_x2=0
        x_previous2=0
        v_previous_y2=0
        y_previous2=0


        v_previous_x3=0
        x_previous3=0
        v_previous_y3=0
        y_previous3=0
    open(12,FILE='result.txt')
        do i=1,5
        v_next_x1=v_previous_x1+const*(abs(x22-x11)/((abs(x22-x11))**2+(abs(y22-y11))**2)&
        +abs(x33-x11)/((abs(x33-x11))**2+(abs(y33-y11))**2)-abs(x11)/((abs(x11))**2+(abs(y11))**2))
        x_next1=x_previous1+tau*v_previous_x1
        v_next_y1=v_previous_y1+const*(abs(y22-y11)/((abs(x22-x11))**2+(abs(y22-y11))**2)&
        +abs(y33-y11)/((abs(x33-x11))**2+(abs(y33-y11))**2)-abs(y11)/((abs(x11))**2+(abs(y11))**2))
                y_next1=y_previous1+tau*v_previous_y1

        v_next_x2=v_previous_x2+const*(abs(x11-x22)/((abs(x11-x22))**2+(abs(y11-y22))**2)&
        +abs(x33-x22)/((abs(x33-x22))**2+(abs(y33-y22))**2)-abs(x22)/((abs(x22))**2+(abs(y22))**2))
        x_next2=x_previous2+tau*v_previous_x2
        v_next_y2=v_previous_y2+const*(abs(y11-y22)/((abs(x11-x22))**2+(abs(y11-y22))**2)&
        +abs(y33-y22)/((abs(x33-x22))**2+(abs(y33-y22))**2)-abs(y22)/((abs(x22))**2+(abs(y22))**2))
                y_next2=y_previous2+tau*v_previous_y2

        v_next_x3=v_previous_x3+const*(abs(x11-x33)/((abs(x11-x33))**2+(abs(y11-y33))**2)&
        +abs(x22-x33)/((abs(x22-x33))**2+(abs(y22-y33))**2)-abs(x33)/((abs(x33))**2+(abs(y33))**2))/me
        x_next3=x_previous3+tau*v_previous_x3
        v_next_y3=v_previous_y3+const*(abs(y11-y33)/((abs(x11-x33))**2+(abs(y11-y33))**2)&
        +abs(y22-y33)/((abs(x22-x33))**2+(abs(y22-y33))**2)-abs(y33)/((abs(x33))**2+(abs(y33))**2))
                y_next3=y_previous3+tau*v_previous_y3

        write(*,*) x_next1,y_next1,v_next_x1,v_next_y1
        !write(12,*) x_next1,y_next1,x_next2,y_next2,x_next3,y_next3


        x11=x_next1
        x22=x_next2
        x33=x_next3
        y11=y_next1
        y22=y_next2
        y33=y_next3

        v_previous_x1=v_next_x1
        x_previous1=x_next1
        v_previous_y1=v_next_y1
        y_previous1=y_next1


        v_previous_x2=v_next_x2
        x_previous2=x_next2
        v_previous_y2=v_next_y2
        y_previous2=y_next2


        v_previous_x3=v_next_x3
        x_previous3=x_next3
        v_previous_y3=v_next_y3
        y_previous3=y_next3
        enddo
    close(12)
    end program
 

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

К тому же во все выражения заряды входят в квадрате. Частично можно было бы сказать, что это скомпенсировано использованием ангстремов в качестве единицы длины, но это $10^{-10}$ м, а не $10^{-10}$ см, у Вас остаются неучтенными два порядка.

Нет, я просто вынес $10^{-10}$ из всех констант, это же никак не повлияет на результат?

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

А теперь внимательно посмотрите на начальные условия. В момент старта электроны №1 и №2 находятся в одной точке. Следовательно, при вычислении их ускорений Вы получите нуль в знаменателе. При этом для вычисления новых координат используются переменные вроде x_previous2, которые исходно все равны нулю.

С начальными условиями я ,действительно, тупанул, но тогда для старта моей системы мне нужны еще скорости электронов, где мне их взять, сколько положить, порядка с? Также мне не понятно какого порядка брать $\tau$?

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

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

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

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

Любую книгу по Фортрану, начиная со стандарта Fortran 90, их и в интернете достаточно много. Идея сводится к тому, что вместо операций с отдельными переменными можно выполнять те же операции с массивами, в итоге не надо будет писать однотипные участки кода для каждой из координат.

. То есть, грамотная запись - это запись каждой переменной в вектор(единичный массив) и далее оперирование с ним, но тогда запись сил усложнится.

[Pphantom]

MathKvant в сообщении #974659 писал(а):

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

Чуть лучше, но часть перечисленных ранее проблем есть и в ней. Например, два комплекта начальных координат, один из которых (правильный) используется для вычисления скоростей, а второй (нулевой) - для изменения координат. В итоге первый шаг как-то срабатывает, а дальше все разваливается.

MathKvant в сообщении #974659 писал(а):

Нет, я просто вынес $10^{-10}$ из всех констант, это же никак не повлияет на результат?

Да, но если Вы сделаете это правильно. Используемые Вами радиусы орбит указаны на самом деле в $\text{см}^{-8}$.

MathKvant в сообщении #974659 писал(а):

С начальными условиями я ,действительно, тупанул, но тогда для старта моей системы мне нужны еще скорости электронов, где мне их взять, сколько положить, порядка с?

Оценить, решив задачу о равномерном движении по окружности вокруг ядра и пренебрегая другими электронами. Для начала это сойдет.

MathKvant в сообщении #974659 писал(а):

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

Тогда давайте чуть детальнее: на кого Вы учитесь? На программиста-быдлокодера или на кого-то другого?

От этого зависит, как относиться к задаче. Если верно первое, то достаточно что-то сдать и забыть (хотя тогда и сама задача, и требование использования Фортрана выглядят несколько странно). Если второе, то Вам лучше забыть о своих привычках.

MathKvant в сообщении #974659 писал(а):

То есть, грамотная запись - это запись каждой переменной в вектор(единичный массив) и далее оперирование с ним

Только не единичный, а одномерный массив.

MathKvant в сообщении #974659 писал(а):

но тогда запись сил усложнится.

Наоборот. В учебнике физики силы покоординатно расписывались или все-таки с помощью векторов?

[MathKvant]

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

Чуть лучше, но часть перечисленных ранее проблем есть и в ней. Например, два комплекта начальных координат, один из которых (правильный) используется для вычисления скоростей, а второй (нулевой) - для изменения координат. В итоге первый шаг как-то срабатывает, а дальше все разваливается.

Спасибо, подправил все, программа считает, но неправильно. У меня траектория каждого электрона эллипс, а меня y(x)- линейная функция, то есть у меня нет практически отрицательных координат. Я $\tau$ подобрал так, чтобы получающиеся скорости электронов совпадали по размерности с начальной скоростью.

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

    program molecul
    implicit none
    real, parameter :: r1=529.17720859! радиус первой орбиты в см*10^(-10)
    real, parameter :: r2=2116.70883436! радиус второй орбиты в см*10^(-10)
    real, parameter :: me=9.10938291e-18! масса электрона в гр*10^(-10)
    real, parameter :: q=4.8032042713! заряд элетрона в сгсэ
    real, parameter :: tau=1e-5! шаг по времени в секундах
    real, parameter :: v_0=2e9! cm/c
    real v_next_x1,v_previous_x1,x_next1,x_previous1,v_next_y1,v_previous_y1,y_next1,y_previous1
    real v_next_x2,v_previous_x2,x_next2,x_previous2,v_next_y2,v_previous_y2,y_next2,y_previous2
    real v_next_x3,v_previous_x3,x_next3,x_previous3,v_next_y3,v_previous_y3,y_next3,y_previous3
    real x11,x22,x33,y11,y22,y33
    real, parameter ::  const=tau*q*q/me
    integer i
    write(*,*) const
    x11=r1! начальные условия
    x22=0
    x33=r2
    y11=0
    y22=r1
    y33=0
        v_previous_x1=v_0
        x_previous1=r1
        v_previous_y1=v_0
        y_previous1=0


        v_previous_x2=v_0
        x_previous2=0
        v_previous_y2=v_0
        y_previous2=r1


        v_previous_x3=v_0
        x_previous3=r2
        v_previous_y3=v_0
        y_previous3=0
    open(12,FILE='result.txt')
        do i=1,1000
        v_next_x1=v_previous_x1+const*(abs(x22-x11)/((abs(x22-x11))**2+(abs(y22-y11))**2)&
        +abs(x33-x11)/((abs(x33-x11))**2+(abs(y33-y11))**2)-abs(x11)/((abs(x11))**2+(abs(y11))**2))

        x_next1=x_previous1+tau*v_previous_x1

        v_next_y1=v_previous_y1+const*(abs(y22-y11)/((abs(x22-x11))**2+(abs(y22-y11))**2)&
        +abs(y33-y11)/((abs(x33-x11))**2+(abs(y33-y11))**2)-abs(y11)/((abs(x11))**2+(abs(y11))**2))

        y_next1=y_previous1+tau*v_previous_y1


        v_next_x2=v_previous_x2+const*(abs(x11-x22)/((abs(x11-x22))**2+(abs(y11-y22))**2)&
        +abs(x33-x22)/((abs(x33-x22))**2+(abs(y33-y22))**2)-abs(x22)/((abs(x22))**2+(abs(y22))**2))

        x_next2=x_previous2+tau*v_previous_x2


        v_next_y2=v_previous_y2+const*(abs(y11-y22)/((abs(x11-x22))**2+(abs(y11-y22))**2)&
        +abs(y33-y22)/((abs(x33-x22))**2+(abs(y33-y22))**2)-abs(y22)/((abs(x22))**2+(abs(y22))**2))

        y_next2=y_previous2+tau*v_previous_y2

        v_next_x3=v_previous_x3+const*(abs(x11-x33)/((abs(x11-x33))**2+(abs(y11-y33))**2)&
        +abs(x22-x33)/((abs(x22-x33))**2+(abs(y22-y33))**2)-abs(x33)/((abs(x33))**2+(abs(y33))**2))

        x_next3=x_previous3+tau*v_previous_x3


        v_next_y3=v_previous_y3+const*(abs(y11-y33)/((abs(x11-x33))**2+(abs(y11-y33))**2)&
        +abs(y22-y33)/((abs(x22-x33))**2+(abs(y22-y33))**2)-abs(y33)/((abs(x33))**2+(abs(y33))**2))

        y_next3=y_previous3+tau*v_previous_y3

        !write(*,*) x_next1,y_next1,v_next_x1,v_next_y1
        write(12,*) x_next1,y_next1,x_next2,y_next2,x_next3,y_next3


        x11=x_next1
        x22=x_next2
        x33=x_next3
        y11=y_next1
        y22=y_next2
        y33=y_next3

        v_previous_x1=v_next_x1
        x_previous1=x_next1
        v_previous_y1=v_next_y1
        y_previous1=y_next1


        v_previous_x2=v_next_x2
        x_previous2=x_next2
        v_previous_y2=v_next_y2
        y_previous2=y_next2


        v_previous_x3=v_next_x3
        x_previous3=x_next3
        v_previous_y3=v_next_y3
        y_previous3=y_next3
        enddo
    close(12)
    end program
 

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

Да, но если Вы сделаете это правильно. Используемые Вами радиусы орбит указаны на самом деле в $\text{см}^{-8}$.

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

Оценить, решив задачу о равномерном движении по окружности вокруг ядра и пренебрегая другими электронами. Для начала это сойдет.

Спасибо, уже учел в программе.

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

Тогда давайте чуть детальнее: на кого Вы учитесь? На программиста-быдлокодера или на кого-то другого?

От этого зависит, как относиться к задаче. Если верно первое, то достаточно что-то сдать и забыть (хотя тогда и сама задача, и требование использования Фортрана выглядят несколько странно). Если второе, то Вам лучше забыть о своих привычках.

Учусь на специалиста в физике конденсированных сред. Хорошо бы еще видеть правильный код, чтобы брать пример с него.

[Pphantom]

MathKvant в сообщении #975420 писал(а):

У меня траектория каждого электрона эллипс,

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

MathKvant в сообщении #975420 писал(а):

а меня y(x)- линейная функция

MathKvant в сообщении #975420 писал(а):

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

Так $\tau$ выражено в секундах или в чем-то другом? Потому что если в секундах, то получаемый результат совершенно естественен - у Вас шаг по времени получается существенно большим, чем характерное время пересечения системы, т.е. после первого же шага взаимодействием электронов с ядром и друг другом можно пренебречь, после чего они летят по прямой.

MathKvant в сообщении #975420 писал(а):

Учусь на специалиста в физике конденсированных сред.

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

Тогда давайте делать все нормально и с самого начала. Запишите уравнения движения (только полноценные - в векторном виде и со всеми физическими коэффициентами, а не тот эрзац, с которого началась эта беседа) и приведите их к виду "ускорение электрона № такой-то равно тому-то". Пока без всяких разностных аппроксимаций, это следующий шаг, до него еще добраться надо.

Параллельно можно найти какую-нибудь книгу по Фортрану (например, Немнюгин C.А., Стесик О.Л. "Современный Фортран", она так себе, но ее очень легко найти в Интернете) и читать.

MathKvant в сообщении #975420 писал(а):

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

Нет уж. Во-первых, так Вы ничему не научитесь, во-вторых, правильный код настолько прост и мал по объему, что "брать пример" с него не удастся - только копировать.

[MathKvant]

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

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

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

Так $\tau$ выражено в секундах или в чем-то другом? Потому что если в секундах, то получаемый результат совершенно естественен - у Вас шаг по времени получается существенно большим, чем характерное время пересечения системы, т.е. после первого же шага взаимодействием электронов с ядром и друг другом можно пренебречь, после чего они летят по прямой.

Да, действительно забыл, что не один электрон, но все равно траектория должна быть замкнутой. Теперь понимаю, что не верно выбрал $\tau$.

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

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

Тогда давайте делать все нормально и с самого начала. Запишите уравнения движения (только полноценные - в векторном виде и со всеми физическими коэффициентами, а не тот эрзац, с которого началась эта беседа) и приведите их к виду "ускорение электрона № такой-то равно тому-то". Пока без всяких разностных аппроксимаций, это следующий шаг, до него еще добраться надо.

Параллельно можно найти какую-нибудь книгу по Фортрану (например, Немнюгин C.А., Стесик О.Л. "Современный Фортран", она так себе, но ее очень легко найти в Интернете) и читать.

Хорошо, опишу свои рассуждения с начала. Пусть есть 3 электрона и покоящееся ядро. На каждый электрон по закону Кулона действует сила притяжения или отталкивания.
$ \\ \sum_{i=1}^3\overrightarrow{F_i}=m_1\cdot \overrightarrow{a_1}\\ k\frac{q_1\cdot q_2}{r_{12}^3\cdot m_1}\overrightarrow{r_{12}}+k\frac{q_1\cdot q_3}{r_{13}^3\cdot m_1}\overrightarrow{r_{13}}+k\frac{q_1\cdot q_0}{r_{10}^3\cdot m_1}\overrightarrow{r_{10}}= \overrightarrow{a_1}\\ k\frac{q_2\cdot q_1}{r_{21}^3\cdot m_2}\overrightarrow{r_{21}}+k\frac{q_2\cdot q_3}{r_{23}^3\cdot m_2}\overrightarrow{r_{23}}+k\frac{q_3\cdot q_0}{r_{30}^3\cdot m_2}\overrightarrow{r_{20}}= \overrightarrow{a_2}\\ k\frac{q_3\cdot q_1}{r_{31}^3\cdot m_3}\overrightarrow{r_{31}}+k\frac{q_3\cdot q_2}{r_{32}^3\cdot m_3}\overrightarrow{r_{32}}+k\frac{q_3\cdot q_0}{r_{30}^3\cdot m_3}\overrightarrow{r_{30}}= \overrightarrow{a_3} $, где 0 соответствует ядру.
Или правильная запись через $\operatorname{grad}$?

[Pphantom]

MathKvant в сообщении #976319 писал(а):

Да, действительно забыл, что не один электрон, но все равно траектория должна быть замкнутой.

Это тоже ниоткуда не следует. Вообще говоря, в этой модели Вы даже финитность движения гарантировать не можете.

MathKvant в сообщении #976319 писал(а):

Или правильная запись через $\operatorname{grad}$?

Да нет, так нормально. Теперь замените $q_i$ и $m_i$ на известные заряды и массы, а затем найдите значение размерного коэффициента в какой-то подходящей системе единиц (выше обсуждались по крайней мере две).

[MathKvant]

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

Теперь замените $q_i$ и $m_i$ на известные заряды и массы, а затем найдите значение размерного коэффициента в какой-то подходящей системе единиц

Если всё перевести в СГСЭ, то заряд $ q$=$4,8\cdot10^{-10}$ и $m=9,1\cdot10^{-28}$, то const=$2,5\cdot10^{8}$СГСЭед/гр$, если еще вынести r, то еще умножить на $10^{10}$
Я нашел ошибки в предыдущей программе:
1) Во-первых, ошибся в производной, когда проекцию сил искал, там $r^{\frac{3}{2}}$.
2)Во-вторых, почему-то в числителях у меня модуль, хотя его быть не должно.
3) Забыл учесть, что заряд у ядра равен сумме всех зарядов электронов

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

    program molecul
    implicit none
    real, parameter :: r1=529.17720859! радиус первой орбиты в см*10^(-10)
    real, parameter :: r2=2116.70883436! радиус второй орбиты в см*10^(-10)
    real, parameter :: me=9.10938291e-18! масса электрона в гр*10^(-10)
    real, parameter :: q=4.8032042713! заряд элетрона в сгсэ
    real, parameter :: tau=1e-6! шаг по времени в секундах
    real, parameter :: v_0=2e9! cm/c
    real v_next_x1,v_previous_x1,x_next1,x_previous1,v_next_y1,v_previous_y1,y_next1,y_previous1
    real v_next_x2,v_previous_x2,x_next2,x_previous2,v_next_y2,v_previous_y2,y_next2,y_previous2
    real v_next_x3,v_previous_x3,x_next3,x_previous3,v_next_y3,v_previous_y3,y_next3,y_previous3
    real x11,x22,x33,y11,y22,y33
    real, parameter ::  const=tau*q*q/me
    integer i
    write(*,*) const
    x11=r1! начальные условия
    x22=0
    x33=r2
    y11=0
    y22=r1
    y33=0
        v_previous_x1=v_0
        x_previous1=r1
        v_previous_y1=v_0
        y_previous1=0


        v_previous_x2=v_0
        x_previous2=0
        v_previous_y2=v_0
        y_previous2=r1


        v_previous_x3=v_0
        x_previous3=r2
        v_previous_y3=v_0
        y_previous3=0
    open(12,FILE='result.txt')
        do i=1,1000
        v_next_x1=v_previous_x1+const*((x22-x11)/(((abs(x22-x11))**2+(abs(y22-y11))**2))**1.5&
        +(x33-x11)/(((abs(x33-x11))**2+(abs(y33-y11))**2))**1.5-3*x11/(((abs(x11))**2+(abs(y11))**2))**1.5)

        x_next1=x_previous1+tau*v_previous_x1

        v_next_y1=v_previous_y1+const*((y22-y11)/(((abs(x22-x11))**2+(abs(y22-y11))**2))**1.5&
        +(y33-y11)/(((abs(x33-x11))**2+(abs(y33-y11))**2)**1.5)-3*y11/(((abs(x11))**2+(abs(y11))**2))**1.5)

        y_next1=y_previous1+tau*v_previous_y1


        v_next_x2=v_previous_x2+const*((x11-x22)/(((abs(x11-x22))**2+(abs(y11-y22))**2))**1.5&
        +(x33-x22)/(((abs(x33-x22))**2+(abs(y33-y22))**2))**1.5-3*x22/(((abs(x22))**2+(abs(y22))**2))**1.5)

        x_next2=x_previous2+tau*v_previous_x2


        v_next_y2=v_previous_y2+const*((y11-y22)/(((abs(x11-x22))**2+(abs(y11-y22))**2))**1.5&
        +(y33-y22)/(((abs(x33-x22))**2+(abs(y33-y22))**2))**1.5-3*y22/(((abs(x22))**2+(abs(y22))**2))**1.5)

        y_next2=y_previous2+tau*v_previous_y2

        v_next_x3=v_previous_x3+const*((x11-x33)/(((abs(x11-x33))**2+(abs(y11-y33))**2))**1.5&
        +(x22-x33)/(((abs(x22-x33))**2+(abs(y22-y33))**2))**1.5-3*x33/(((abs(x33))**2+(abs(y33))**2))**1.5)

        x_next3=x_previous3+tau*v_previous_x3


        v_next_y3=v_previous_y3+const*((y11-y33)/(((abs(x11-x33))**2+(abs(y11-y33))**2))**1.5&
        +(y22-y33)/(((abs(x22-x33))**2+(abs(y22-y33))**2))**1.5-3*y33/(((abs(x33))**2+(abs(y33))**2))**1.5)

        y_next3=y_previous3+tau*v_previous_y3

        !write(*,*) x_next1,y_next1,v_next_x1,v_next_y1
        write(12,*) x_next1,y_next1,x_next2,y_next2,x_next3,y_next3


        x11=x_next1
        x22=x_next2
        x33=x_next3
        y11=y_next1
        y22=y_next2
        y33=y_next3

        v_previous_x1=v_next_x1
        x_previous1=x_next1
        v_previous_y1=v_next_y1
        y_previous1=y_next1


        v_previous_x2=v_next_x2
        x_previous2=x_next2
        v_previous_y2=v_next_y2
        y_previous2=y_next2


        v_previous_x3=v_next_x3
        x_previous3=x_next3
        v_previous_y3=v_next_y3
        y_previous3=y_next3
        enddo
    close(12)
    end program
 

[Pphantom]

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