Ещё один путь к Земле

Представляем вашему вниманию оригинальную реализацию программ серии "Путь к земле" в виде веб-страниц со скриптами на сайте Станислава Володарского.

"На этом сайте есть всё необходимое, чтобы повторить перелёт «Кон-Тики» с Луны на Землю и решить многие другие задачи космонавигации."

Кроме программ-симуляторов, лунолетчикам может потребоваться книга "Задачи для юного космонавта".

Комментарии

Доводы против EQU понятны.

Возможность авторазмещения регистров рассматривалась, хотя пока не реализована. Особых сложностей ввести её нет.

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

Кстати, mk.exe работает под wine. Нужно только не забыть последовательный порт в /dosdevices прописать.

Отладка на уровне MKL файла. Тогда и EQU и символические метки засияют с новой силой.

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

Насколько я понимаю, логика вызова Дорманд-Принса будет отличаться. Вместо полного и двух половинных расчётов мы будем считать f(p0, t, ∆t) дважды — аргументы одинаковы, а вот формулы теперь будут разные. При этом на повторном расчёте можно будет использовать предыдущие pN' . Вот только коэффициенты будут во взвешенной сумме производных будут далеко ее так хороши.

В-общем, мы подлетаем к вашей вершине.

Коэффициент для РК будет 16/15.

Общая формула: coeff = 2^k/(2^k - 1), где k - порядок метода (4 для РК, 1 для Э)

Перерасчёт dt при ошибке:

new_dt = 0.9 * dt / (rer^(1/k))

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

Коэффициент для методов высшего порядка можно смело приравнять единице. Множитель 0.9 его эффективно съедает (особенно после извлечени корня порядка k).

Добавил перерасчёт коэффициента и ∆t, в зависимости от метода. Не проверял. У хостинга какие-то проблемы, выложил на народ:
lun-103.zip.html

Смелость, думаю, излишня, т.к. коэффициент также используется для оценки ∆t, а большинство недостатков команды Fxy исправлено. Мы и так уже схитрили с методом, введя dt_min. Успокаивает лишь надежда, что наша система действительно так хороша, как вы её описали и обычно прощает эту хитрость.

Я посмотрел вариант Рунге-Кутты. Он не заработает. Нет скоростей в переменных p1' - p4'. Там сейчас только ускорения. Нужны также и скорости.

// размерности векторов не совпадают
x, y, vx, vy := mp + dt * (ax, ay)

// а так совпадают
x, y, vx, vy := mp + dt * (vx, vy, ax, ay)

Предлагаю разбираться в коде, а не придираться к комментариям. Программа у меня же считала. У вас она работает?
C размерностями всё в порядке точно.

Возможно, что комментарий неудачен. Я так кратко сформулировал вот такое вычисление, которое осуществляется подпрограммой PlusAnaDT:
vx := mp.vx + ax * dt
vy := mp.vy + ay * dt
x := mp.x + mp.vx / mp.y * dt
y := mp.y + mp.vy * dt

Теперь о переменных p1'-p4' — это производные? Я могу отвести для них ещё регистры, но откуда в них возьмётся скорость? Дело в том, что вы предложили не передавать аргумент dt блоку "ФИЗИКА", то есть g(). Запомнить скорости в p1'-p4' невозможно — функция g() их не в состоянии рассчитать и вернуть, не зная dt.

Расчёт скоростей происходит при вычислении первого аргумента функции g(). Занимается этим подпрограмма PlusAnaDT, уже после вычисления ускорений. Формулы см. выше. Для этого я ей передаю dt, которое известно лишь функции f().

Всё в соответствии с вашим псевдокодом, разве что я осмелился ввести ещё одну подпрограмму PlusAnaDT для упрощения вычисления первого параметра g().

Подсчёт скоростей был вынесен в конец блоков Эйлера и РК, на основе взвешенного ускорения (производных, в вашей терминологии). Если скорости нужно подсчитывать в блоке "ФИЗИКА", то прошу прописать явно, какое dt этому блоку, то есть функции g() в каждом из четырёх вызовов передавать.

Если вы всё ещё видите ошибку, которая отсутствует в псевдокоде, пожалуйста, укажите на неё точнее.

Программа у вас и у меня считала и сейчас считает великолепно (я её тестировал-перетестировал, всё в порядке, а как она летает по эллиптическим орбитам!).

Когда используется метод Эйлера скорости берутся из начала шага. Поэтому при экстраполяции можно воспользоваться исходными скоростями. Но РК делает экстраполяции не только из начала шага, но и из середины и из конца. Там скорости должны быть другими - проэкстраполированными. Следовательно, я предлагаю, чтобы блок "физика" возвращал не два значения (ускорения), а четыре (ускорения и копии скоростей, которые он получил как аргументы).
Блок PlusAnaDT должен получать четыре параметра - сейчас он две скорости всегда берёт из начала шага.
Это отбрасывает нас назад, к Эйлеру. То есть, считать будет, но с точностью первого порядка, так как скорости считаются постоянными.

Переменные p1' - p4' действительно производные. Это производные координат и скоростей. Производные скоростей - ускорения, производные координат - скорости и их тоже надо учитывать.

Таким образом g должен возвращать ускорения (они сейчас считаются правильно) и копии скоростей. g по прежнему не получит dt.

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

> g должен возвращать ускорения (они сейчас считаются правильно) и копии скоростей

У меня от такого глаза полезли на лоб. Ладно, сделаю. Раз хотите, в чём вопрос. :-)

Но вам придётся расписать покомпонентно следующие две стенограммы (первая для PlusAnaDT, вторая для return):

x, y, vx, vy := mp + dt * (vx, vy, ax, ay)

x_new, y_new, vx_new, vy_new := mp + dt * (vx, vy, ax, ay)

экстраполяция выглядит так:
new.x = x + dt * p'.vx
new.y = y + dt * p'.vy
new.vx = vx + dt * p'.ax
new.vy = vy + dt * p'.ay

Индекс после p я нарочно опустил.
Если задуматься, то окажется, что метод Эйлера от такой переделки не поменяется, так как в нём p'.vx == vx, p'.vy == vy.
И первый прогноз РК тоже не изменится. А вот второй получит изменившиеся скорости и учтёт их изменения в расчётах.

Это не я этого хочу, это мы решаем диффур второго порядка (координаты, скорости, производные).
Вас ведь не удивило, что на взлёте после первого шага метода Эйлера, вертикальная скорость положительна, а высота нулевая (это я про баг с посадкой на взлёте)! И только второй шаг сдвигает лунолёт с места. То есть, метод Эйлера переносит влияние ускорения на скорость за один шаг, а на координату за два. Зато он прост как доска, и работает. И РК мы заставим работать, и в нём ускорение будет влиять на все величины сразу.

А ничего, что координата x у нас — в радианах? И точно ли, что скорость (vx, vy) не влияет на (new.x, new.y) — только (p'.vx, p'.vy). Кстати, раз уж мы перешли на ваши обозначения, какой псевдокод вы расшифровываете — вот этот?

(new.x, new.y, new.vx, new.vy) = (x, y, vx, vy) + dt * p'

Про Эйлера можно забыть, до этапа оптимизации. Я его уже отделил от Рунге-Кутта, чтобы в связи с новыми веяниями не покорёжить.

Хостинг, кстати, починили. Обновил файлы и там — сейчас там та же версия, что и на народе.

Для удобства сведу весь псевдокод вместе. Проверьте, пожалуйста. Я постарался устранить все ошибки и опечатки. Вы были правы, скорости не совсем копируются.

Итак мы решаем систему диффуров. Вот она в векторном виде:

p' = g(p, t)

здесь p = (x, y, vx, vy),
p' = (x', y', vx', vy') = (vx, vy, ax, ay)

Вот она в виде системы:

x' = vx / y
y' = vy
vx' = rax(t) - vy * vx / y
vy' = ray(t) + vx * vx / y - gc / y^2

Вот она в виде функции:

g(x, y, vx, vy, t) {
    x' = vx / y
    y' = vy
    vx' = rax(t) - vy * vx / y
    vy' = ray(t) + vx * vx / y - gc / y^2
    return (x', y', vx', vy')
}

Функция h выполняет линейную экстраполяцию:

h(x, y, vx, vy, dt, x', y', vx', vy') {
    new_x  = x  + dt * x'
    new_y  = y  + dt * y'
    new_vx = vx + dt * vx'
    new_vy = vy + dt * vy'
    return (new_x, new_y, new_vx, new_vy)
}

g и h в векторном виде:

g(p, t) {
    p'.x' = p.vx / p.y
    p'.y' = p.vy
    p'.vx' = rax(t) - p.vy * p.vx / p.y
    p'.vy' = ray(t) + p.vx * p.vx / p.y - gc / p.y^2
    return p'
}

h(p, dt, p') {
    new_p = p + dt * p' // операции над векторами
    return new_p
}

Так мы их объединяем:

f(p, t, dt) {
    p' = g(p, t)
    new_p = h(p, dt, p')
    return new_p
}

Опишем шаг метода Эйлера (функция f его реализует полностью):

euler(p, t, dt) {
    return f(p, t, dt)
}

Или так:

euler(p, t, dt) {
    return h(p, dt, g(p, t))
}

Опишем шаг метода Рунге-Кутты:

rungeKutta(p, t, dt) {
    p1' = g(p              , t       ) // начало
    p2' = g(h(p, dt/2, p1'), t + dt/2) // середина
    p3' = g(h(p, dt/2, p2'), t + dt/2) // ещё раз
    p4' = g(h(p, dt  , p3'), t + dt  ) // конец шага!

    return h(p, dt, (p1' + 2 * (p2' + p3') + p4') / 6);
}

Делаем шаг и вычисляем ошибку: (m - функция euler или rungeKutta):

stepAndError(p, t, dt) {
    p1 = m(p, t, dt)
    p21 = m(p, t, dt/2)
    p2 = m(p21, t + dt/2, dt/2)
    return (p1, |p1 - p2|)
}

Теперь всё довольно чётко расписано. Если x' у нас радианы в секунду, действительно можно делать линейную экстраполяцию с координатой x. Вообще, спасибо за функцию h(). Приятно, когда математику даёт человек, умеющий программировать.

Единственное, что меня настораживает и настораживало раньше — g(p, t) не использует координату p.x . Но, видимо, это особенности нашей физической модели. Также я не совсем понимаю знак минус у Кориолиса. Но это вряд ли ваши сложности. Понять Кориолиса мне удалось лишь один раз в жизни, потом опять забылось. Это вообще одна из моих личных целей всего проекта. Центробежка плюс, гравитация минус — здесь всё очевидно. Отдельные ошибки по скорости и координатам я распишу сам, этот код уже есть.

День у меня занят, но ночью буду кодировать для ЭКВМ. Перейду на эти обозначения. Многие регистры теперь будут распределены по-новому. Но Дорманд-Принс, насколько я понял, так легче будет встроить.

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

Пусть vx0 > 0 и vy0 == 0, свободный полёт, даже гравитации нет. Отправляем корабль по касательной к окружности. При этом vy начинает расти. Модуль скорости должен сохранится (тяги нет, тяготения тоже). Следовательно vx^2 + vy^2 = const. По мере удаления от тела вертикальная скорость растёт, следовательно горизонтальная должна уменьшаться (с течением времени vy подберётся почти к vx0, а vx станет почти нулём). Обратите внимание, что физики никакой тут нет - только геометрия для движения с постоянной скоростью в прямоугольной системе координат.

Какие причины заставляют расти vy? Центробежное ускорение (фиктивное). Нужен аналогичный механизм для уменьшения vx. Это кориолисово ускорение (тоже фиктивное).

Если у вас vx != 0, то центробежное ускорение положительно. Пусть vx > 0 и vy > 0. Под действием центробежного ускорения vy растёт, из сохранения полной скорости следует, что кориолисово ускорение должно уменьшать vx.
Отсюда правило знаков: vx = vx - dt * (vx * vy / y). Кориолисово ускорение всегда с минусом.

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

Для ДП придётся усложнить функцию h. Но сперва РК.


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

Эх, не хотелось затевать это до написания «Лунолёта». Ещё не встречал физика, который бы это понимал — ни в физматшколе, ни на физфаке МГУ. Слишком мало их жило (или вживалось) в неинерциальные системы отсчёта, большинство считает эти ускорения математической абстракцией. Типа математики такой трюк выдумали, физического смысла не сказали — значит его как бы и нет. :-)

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

Твоему же телу и всем предметам на борту хочется продолжать двигаться по прямой инерции (которая субъективно приводит ко всё большему увеличению высоты над поверхностью небесного тела, будто нас всех отталкивает от небесного тела какая-то сила, не зависящая от массы) и пилот будет ощущать эту центробежку v²/r из-за насильного загибания траектории до окружности h = const. Прижмёт к потолку кабины.

Теперь учтём гравитацию, она будет оттягивать нас (мы внутри кабины корабля, на корпус которого действует реактивная сила) вниз, к небесному телу. В идеальном случае гравитация уравновесит центробежку и наступит невесомость. Никакая третья сила, навроде реактивной тяги, для реализации такой неинерциальной системы отсчёта уже не будет нужна — мы вышли на орбиту и движемся с первой космической.

С помощью постоянной тяги мы можем двигаться по орбите h=const с любой другой скоростью — больше или меньше первой космической. В этом случае мы будем ощущать в кабине искусственную гравитацию (точнее, искусственную силу тяжести = сила гравитации + силы инерции) — либо вниз, вплоть до полного ускорения свободного падения при нулевой орбитальной скорости. Либо отрицательную вверх, если пойдём крутиться со скоростью больше первой космической.

Вот также расписать Кориолиса, на субъективных ощущениях в желудке космонавта, я пока не могу. :( Однажды я эту силу инерции «поймал», дотошно разбирая физический смысл каждой из формул при выводе формул преобразования координат, но для закрепления нужен опыт полётов.

> Твоему же телу хочется продолжать двигаться по прямой ... и оно будет ощущать эту центробежку v²/r из-за насильного загибания траектории до окружности h = const. Прижмёт к потолку кабины.

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

> ... (точнее, искусственную силу тяжести = сила гравитации + силы инерции) ...

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

В школе нам объясняют, что есть сила тяготения и она направлена вниз. Но с точки зрения наших ощущений это не мы давим на землю, это она на нас давит. И она действительно на нас давит. Это сила реакции опоры и она направлена вверх. Но мы не улетаем, значит, говорят физики, сумма сил равна нулю. Тогда мы вводим силу тяготения, направленную вниз - физики довольны. С течением времени у нас возникает впечатление, что мы ощущаем силу тяготения и она направлена вниз. Обычный человек забывает про силу реакции, но твёрдо помнит про гравитацию.

Резюме. Что-то прижимает нас к земле и мы думаем, что чувствуем тяготение. Но мы не чувствуем тяготение.

И так, у нас есть привычка заменять реальную силу на противоположную ей, но неощутимую на практике. Перебираемся в космос. Корабль движется по круговой орбите на скорости больше первой космической. На нас давит стенка корабля. Что это? Мы же в относительном покое. Привычная со школы схема, требует чтобы существовала некоторая сила, которая прижимает нас к стенке (как тяготение прижимает нас к земле). Мы придумываем такую силу, называем её центробежной и заблуждение готово.

Центробежная сила не может вас ни к чему прижать (как и гравитация в космосе). Если вы в космосе чувствуете тяжесть, вы просто ускоряетесь в противоположном направлении под действием силя реакции опоры.

Схема проста - сила реакции опоры и ускорение. Но мы же покоимся на опоре! Покоимся! Какое ускорение! Покой невозможен без равновесия сил!
Для объяснения покоя сознание требует вторую силу - аналог гравитации. Это привычно, но не правильно.

Вот примерно так встреченные мною физики и рассуждали. Дескать, физический смысл имеют лишь инерциальные системы отсчёта. А все остальные — «неправильные» (грешные, порочные, недостойные настоящих мужчин и т.д.) :-))

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

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

Силы инерции вполне реальны для человека, который твёрдо решил связать свою систему отсчёта с кораблём, то есть считать его неподвижным. От них можно умереть, ими можно зачать ребёнка. :-) И да, если тело покоится, силы инерции действительно уравновешиваются реакцией опоры.

В инерциальной системе отсчёта эти же явления объясняются по другому, не прибегая к силам инерции.

>Корабль движется по круговой орбите на скорости больше первой космической. На нас давит стенка корабля. Что это?

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

Так трудно двигаться по орбите со скростью, большей первой космической (для данной высоты) без включённых двигателей. Будет сносить от планеты в космос. Двигатели же «прижимают» корабль к орбите, увеличивая центростремительное ускорение и не давая ему улетать, не взирая на завышенную скорость.

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

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

Давайте разберёмся с терминологией.

Я имел в виду движение по строго круговым орбитам со скоростью V, превышающей √(K/R). Где K — гравитационная постоянная планеты, а R — радиус орбиты. При этом V может спокойно превышать и вторую, и даже третью космическую. :-)

Чтобы выполнить такой трюк на современных реактивных кораблях, двигатель включить придётся.

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

И не только на борту — управляя силами инерции, наш реактивный двигатель может притянуть или отфутболить ничего не подозревающую Землю так, что мало не покажется! :-))) Вот только такое же интуитивное представление о Кориолисе, с точки зрения пилота корабля, требует разработки. :(

Выдалось время поработать.

Некоторые новые имена переменных сохранить не удалось. Дело в том, что в архитектуре МК-161 все переменные глобальные и у них должны быть разные имена, точнее регистровые номера. С другой стороны, чтобы не плодить их до бесконечности, некоторые удалось объединить — досрочная оптимизация, хотя я это и не приветствую. Например, h() пишет результат в q, который является аргументом g().

Если эта предварительная оптимизация будет мешать при ДП, разнести регистры — не проблема.

Вот, что получилось — проверяйте. Код на ЯМК пока не выкладываю, т.к. ещё осталось закодировать Рунге-Кутта. Он будет сложнее, хотя и очевидней. :-) Надеюсь ночью закончить.

g(q, tt) {
    p'.x' = q.vx / q.y
    p'.y' = q.vy
    p'.vx' = rax(tt) - q.vy * q.vx / q.y
    p'.vy' = ray(tt) + q.vx * q.vx / q.y - gc / q.y^2
    return p'
}


h(mp, dt, p') {
    q = mp + dt * p' // операции над векторами
    return q
}


euler(mp, mt, dt) {
    return h(mp, dt, g(mp, mt))
}

rungeKutta(mp, mt, dt) {
    p1' = g(mp              , mt       ) // начало
    p2' = g(h(mp, dt/2, p1'), mt + dt/2) // середина
    p3' = g(h(mp, dt/2, p2'), mt + dt/2) // ещё раз
    p4' = g(h(mp, dt  , p3'), mt + dt  ) // конец шага!

    return h(mp, dt, (p1' + 2 * (p2' + p3') + p4') / 6);
}

if (iMetod == 0) m = euler
if (iMetod == 1) m = rungeKutta

stepAndError(p0, t, ∆t) { 
    p1 = m(p0, t, ∆t)
    p21 = m(p0, t, ∆t/2)
    p2 = m(p21, t + ∆t/2, ∆t/2)
    return (p1, |p1 - p2|)
}

Вот предварительная сводка изменений в регистрах, скорее для комплекта. Для понимания нового псевдокода она не требуется.

Переименовать
p0 (R26…R29 = x, y, vx, vy)
p1 (R36…R39 = xy, v = x, y, vx, vy)
pp1 → p21 (R46…49 = x, y, vx, vy)
pp2 → p2 (R50…53 = xy, v = x, y, vx, vy)
x, y, vx, vy → q (R30…R33 = x, y, vx, vy)
p’ (R88,89,78,79 = x’, y’, vx’, vy’)
p1’ (R80,81,70,71 = x’, y’, vx’, vy’)
p2’ (R82,83,72,73 = x’, y’, vx’, vy’)
p3’ (R84,85,74,75 = x’, y’, vx’, vy’)
p4’ (R86,87,76,77 = x’, y’, vx’, vy’)
mp (R20…23 = x, y, vx, vy)
tt (R66)
mt (R56)

Ошибок не заметил.

Хочу чтобы вы подумали над передачей параметров через "указатели". Функция может получать в стеке (или в инкрементных регистрах) номера регистров из которых она будет читать или в которые она будет писать:

26 П4
80 П5
46 П6
dt П7
ПП :h

:h
КИП4 КИП5 ИП7 * + КП6
КИП4 КИП5 ИП7 * + КП6
КИП4 КИП5 ИП7 * + КП6
КИП4 КИП5 ИП7 * + КП6
В/О

В ДП нам без таких "указателей" туго придётся.

Хорошо. Продолжаю кодировать Рунге-Кутта по новому псевдокоду. Кстати, уже заметил ошибку: в методе Рунге-Кутта параметр dt в h() не всегда равен параметру dt в m(). Иногда он равен dt/2, поэтому придётся найти другое обозначение (и регистр) этому аргументу функции h().

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

Таблицы удобнее располагать в десятичных регистрах старше 99, например 126, 180, 246 — младшие регистры мы уже почти все использовали. Вызов подпрограмм в нашей системе кодирования а-ля Huge будет осуществляться вне зависимости от страницы, трёхбайтовая команда ПП с префиксом Р. В традиционной мнемонике ваш код будет выглядеть примерно так:

     126 П4
     180 П5
     246 П6
     dt П7
     РПП h
     …
h:
     В↑
     3 П0
Loop:
     КИП4 КИП5 ИП7 * + КП6
     PFL0 Loop
     В/О

Выложил честного Рунге-Кутта по тому же адресу:
http://the-hacker.ru/2011/lun-103.mkl
http://the-hacker.ru/2011/lun-103.mkp
http://the-hacker.ru/2011/lun-103.htm

Не тестировал. Избавился от хэка с таблицей адресов. В распределении первой сотни регистров осталось несколько «дыр», которые можно использовать для Дорманд-Принса. Среди них — все регистры R0…R9 (удобные для косвенной адресации) и R90…R99 (в принципе, тоже можно использовать для косвенной адресации).

Осталось более 900 десятичных регистров и более 9000 шагов. Вполне достаточно для девяти таких же программ. :-)

Спасибо. Летать буду уже завтра.

Спокойной ночи! Залил новую версию 0.8

Рунге-Кутт почему-то зацикливается. Эйлер же думает долго, но даёт правильный ответ.
Где-то кто-то из нас намудрил. :(

Баг:

		PRM88 PM82 PRM89 PM83 PRM78 PM72 PRM79 PM23	; p2' := p'

Должно быть

		PRM88 PM82 PRM89 PM83 PRM78 PM72 PRM79 PM73	; p2' := p'

Спасибо. Исправил и закачал обратно. Теперь не зависает, но у меня результаты при R57=0

:45 r/s 90 r/s 5 r/s
 Г                0097
-----------------------
T  197.58367356798     
Z  155.60153254904     
Y  1677.6049183167     
X  742.13036098214     
 Угол тяги, градусы ? 

и R57=1

:45 r/s 90 r/s 5 r/s
 Г                0097
-----------------------
T  241.29668611487     
Z  52.582354068047     
Y  1677.5525481987     
X  850.69274496867     
 Угол тяги, градусы ? 

различаются. :(

Так что, проверяя новую версию, не удаляйте свою работающую.

Возможно, у вас Эйлер упирается в dt_min. Давайте уберём dt_min и сделаем аварийный останов при нулевом dt.
И не забудьте поднять точность.

Это я глупость написал. Эйлер у вас на правду похож, а РК нет.

Эйлер как раз даёт числа, совпадающие с сайтом.

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

Ваш результат для РК возникает при задании команды

1 90 5

Ищите ошибку в управлении.

Спасибо, разобрался. Действительно ошибка оператора — не смог правильно набить команду PM57. Получилось PPM 57 и число 45 воспринялось, как номер регистра R5745. Команда на двигатель действительно ушла с углом 1 градус (код РК), а расчёт произошёл по Эйлеру.

Хорошо, у меня логи ведутся. А вы прямо Шерлок Холмс. :-)

Выложил версию 0.10 — так, в поисках несуществующей ошибки чуть причесал исходный код. Работает нормально.

Стоит ли менять ошибки, зашитые в блоке инициализации? По умолчанию же стоит Эйлер, а он медленный. Опять же, на железке будут доступны всего 12-14 десятичных разрядов.

Да, стоит поменять и ошибки и метод. Практически Эйлер бесполезен. Его надо оставить для самоконтроля на небольших манёврах. Но летать надо уже на РК. И ДП будем сравнивать уже с РК.

Сейчас у нас они вот так задаются:

                0,001 PM34 PM35         ; exy, ev
                0,000001 PM43           ; dt_min

                Cx PM57                 ; iMetod := 0 (Эйлер)

На что менять? Сразу же проверю на МК-161, т.к. там не должно тормозить.

Метод надо заменить на РК. Думаю, что можно ошибку взять от 10^-6 до 10^-9.

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

1. Берут манёвр, вроде нашего тестового 45 90 5. А лучше сжесть три с половиной тонны сжечь секунд за 100.

2. Считают манёвр (целиком манёвр, не один шаг) два раза с оценками ошибками (e), отличающимися в два раза. Считают ошибку - разницу между векторами в конце манёвра (можно на глазок прикинуть покомпонентно).

Получается функция которая по оценке ошибки получает саму ошибку для конкретного манёвра.

s(e) = |maneur(e) - maneur(e/2)|
- здесь maneur - возвращает вектор скоростей и координат после тестового манёвра (e - оценка ошибки),
|| - стандартный расчёт ошибки для двух векторов.

3. Функцию s считают для начального e (0,001), потом для в два раза меньшего, в четыре, в восемь и т.п.

4. Пока с уменьшением e функция s ведёт себя линейно, всё в порядке (должно быть примерно s(e/2) = s(e)/2).

5. Начиная с некоторого значения ee s(e) перестанет уменьшаться. Мы упёрлись в "шум", наведённый округлениями.

6. Надо подобрать ближайший больший десятичный порядок. Увеличить его на один порядок. Это будет наш нижний предел для e.

7. Если машинка с таким e работает достаточно быстро (что такое быстро решаете вы), то его можно оставить.

8. Если калькулятор тормозит, то e можно увеличить (надо выбирать между точностью и комфортом игрока - это тоже решаете вы).

Важно понимать что в нашей модели время расчётов (число шагов) пропорционально времени манёвра (при неизменных ускорениях) и пропорционально корню четвёртой степени из ускорений (при неизменном времени манёвра):

<число шагов> = const * DT * sqrt(sqrt(a + ga)) - второй множитель взят сною с потолка :)

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

Понятно. Ну, это может сделать любой затившийся наблюдатель с машинкой. Когда мы Дорманд-Принс доделаем.

Пока сделал Рунге-Кутта методом по умолчанию и поставил 10–6, половину 12-разрядной сетки, а там посмотрим. У меня даже в этом значении сомнения, т.к. радиус Луны много отъедает точности от координаты y, а ведь по ней рассчитывается гравитация и центробежка.

Выложил версию 0.11, на том же месте. Заодно задал более точные значения радиуса и гравитационной постоянной Луны. Глупо претендовать на шесть знаков, опираясь на наивные 1,62 м/с².

Стабильный Рунге-Кутт (версия 0.11) можно скачать отсюда. Вряд ли мы его менять будем, по плану дальше пойдёт работа над Дорманд-Принсем:
http://the-hacker.ru/2011/lun-103.mkl

После первого С/П исследуемую ошибку e (сейчас там 1 ВП –06) вводить в регистры R34 и R35. Никто не мешает её ввод прописать в коде, конечно.

Вектор скоростей брать из R28 и R29. Полярные координаты из R26 и R27 — их разница считается, как √((x1*y1 - x2*y2)² + (y1 - y2)²). Это единственная неочевидная формула. Если не интересно считать ручками, подсчёт погрешности s(e) = |maneur(e) - maneur(e/2)| можно запрограммировать, скажем, с адреса 8000, используя регистры R900…R999 для запоминания результатов прошлого манёвра.

Это по программе. Метод же подсчёта оптимальной точности e написал выше Станислав, снабдивший нас использованной математикой. Не уверен, что всё автоматизируется. Скорее всего, для подбора нужной точности придётся сделать 20-40 прогонов программы и вручную проанализировать, где разрядная сетка начнёт давать сбой. Для Дорманд-Принса придётся повторить, но уже будет опыт, на который можно опираться.

Будет забавно, если из-за ограниченности разрядной сетки Рунге-Кутт окажется точнее. :-)

Нужные регистры, как понимаю, инициализируются автоматически, а тестовый манёвр взять {45, 90, 5} - для сравнения с вашими данными.

Если можно, чуть подробнее про алгоритм вычисления maneur(e) - с конкретными номерами регистров.

это эвфемизм для содержимого четырёх регистров R26, R27, R28 и R29 после команды манёвра 45 С/П 90 С/П 5 С/П (возможно Станислав подберёт команду и начальные условия получше), отработанной с точностью e, заданной в регистрах R34, R35.

Как писал AtH:

После манёвра скорость будет в R28 и R29. Полярные координаты в R26 и R27.
Летаем два раза с разными оценками ошибок e и e/2. Получаем два вектора вида (vx, vy, x, y).

Ошибки считаем так: ev = sqrt((v1x - v2x)^2 + (v1y - v1y)^2)
exy = sqrt((x1*y1 - x2*y2)^2 + (y1 - y2)^2).

Пара этих чисел - наша ошибка (результат функции maneur).

Повторяем процесс для e' = e/2 и e'/2.

Снова считаем ошибки ev' и exy'. Верно ли, что ev' = ev / 2 и exy' = exy / 2 (примерно, с точностью до 10%)?

Если да, то исходная оценка ошибки хорошая. Идём дальше вниз (делим оценку ошибки на 10). И так отыскиваем самую маленькую хорошую оценку ошибки.

ev и exy не обязательно рассчитывать. Достаточно оценивать покомпонентные разницы x1 - x2, y1 - y2, vx1 - vx2, vy1 - vy2. Пока они убывают линейно и ev с exy будут убывать линейно. Когда хотя бы одна из них перестала убывать линейно, минимальная хорошая оценка ошибки найдена.

Версия 0.11 с уточнёнными параметрами Луны, выложенная на http://the-hacker.ru/2011/lun-103.mkp работает на МК-161. К сожалению, основной экран версией 1.21 MK.EXE считать не удалось ("Скорость 9600 бит/с||Строка экрана не считана||Ошибка при считывании"), поэтому переписываю данные вручную (манёвр 45/90/5):

T 197,562636769
Z 155,601546781
Y 1679,4968495348
X 742,0808

Время работы 16 секунд. Ниже тоже самое на эмуляторе.

:45 r/s 90 r/s 5 r/s
 Г                0110
-----------------------
T  197.56263676293     
Z  155.60154677624     
Y  1679.4968495332     
X  742.08080363477     
 Угол тяги, градусы ?

Можно заметить, как большой радиус Луны обрезал на МК-161 точность высоты. От скоростей тоже осталось лишь 12 знаков, т.к. они взяты из регистров. Конечно, это лучше, чем на ПМК. Там высота вообще округлялась до десятых.

А вот тот же расчёт на кастрированном ВК-6 (точность записи числа в десятичных регистрах снижена до 12 знаков):

:45 r/s 90 r/s 5 r/s
 Г                0110
-----------------------
T  197.562636768       
Z  155.60154678        
Y  1679.4968495348     
X  742.0808            
 Угол тяги, градусы ?   

Видно, что такой вариант ВК-6 достаточно неплохо передаёт точность МК-161. Максимальная разница — одна единица младшего разряда. Добиться большего сходства можно разве что скопировав один-в-один bcd-реализацию арифметики и функций.

Отличный результат!

Как меняется время прогона при изменении точности в 10 раз в обе стороны?

1 ВП –05 РП34 РП35 даёт 7,3 секунды — почти реальное время.

T 197,562636912
Z 155,601546908
Y 1679,4968495396
X 742,08079

1 ВП –07 РП34 РП35 даёт 16,1 секунды.

T 197,562636751
Z 155,601546765
Y 1679,4968495348
X 742,0808

Время не очень точное, я вручную засекал. Как насчёт того, чтобы перейти к Дорманд-Принсу? Если помните, именно он меня так заинтересовал. Все эти Эйлеры и Рунге-Кутты были сделаны исключительно по вашему настоянию — мне было проще написать код, чем возражения. :-)

Удивительно что времена счёта для 10^-6 и 10^-7 одинаковы. Мы "вышли" из разрядной сетки калькулятора?

ДП скоро будет. На Э и ДП я настаивал, так как боялся сложных багов на стыке метода, "физики" и управления. Сейчас у нас всё отлажено кроме метода.

План работ:
1. Объяснение работы табличного метода.
2. Изменения в верхнем уровне программы.

1.
Вот наш Рунге-Кутта:

g(q, tt) {
    ...
    return p'
}

h(mp, dt, p') {
    return mp + dt * p' // операции над векторами
}

rungeKutta(mp, mt, dt) {
    p1' = g(mp              , mt       ) // начало
    p2' = g(h(mp, dt/2, p1'), mt + dt/2) // середина
    p3' = g(h(mp, dt/2, p2'), mt + dt/2) // ещё раз
    p4' = g(h(mp, dt  , p3'), mt + dt  ) // конец шага!
    return h(mp, dt, (p1' + 2 * (p2' + p3') + p4') / 6);
}

Вот таблица для РК один раз с конкретными числами, второй - в общем виде:

0   |                     0  |
1/2 | 1/2                 c2 | c21
1/2 | 0   1/2             c3 | c31 c32
1   | 0   0   1           c4 | c41 c42 c43
----------------------   ----------------------
    | 1/6 1/3 1/3 1/6        | c51 c52 c53 c54

Строка таблицы определяет как по предыдущим производным построить новую.
Нам понадобится функции h2 - h5:

h2(mp, dt, c21, p1') {
    return mp + dt * (c21 * p1')
}
h3(mp, dt, c31, p1', c32, p2') {
    return mp + dt * (c31 * p1' + c32 * p2')
}
h4(mp, dt, c41, p1', c42, p2', c43, p3') {
    return mp + dt * (c41 * p1' + c42 * p2' + c43 * p3')
}
h5(mp, dt, c51, p1', c52, p2', c53, p3', c54, p4') {
    return mp + dt * (c51 * p1' + c52 * p2' + c53 * p3' + c54 * p4')
}

Метод РК в табличном виде:

rungeKutta(mp, mt, dt) {
    p1' = g(   mp                                   , mt          ) // начало
    p2' = g(h2(mp, dt, c21, p1'                    ), mt + c2 * dt) // середина
    p3' = g(h3(mp, dt, c31, p1', c32, p2'          ), mt + c3 * dt) // ещё раз
    p4' = g(h4(mp, dt, c41, p1', c42, p2', c43, p3'), mt + c4 * dt) // конец шага!
    return h5(mp, dt, c51, p1', c52, p2', c53, p3', c54, p4');
}

Продолжение следует...

Таблица для Дорманда-Принса:

0    |
1/5  |     1/5
3/10 |     3/40          9/40
4/5  |    44/45        −56/15       32/9
8/ 9 | 19372/6561   −25360/2187  64448/6561   −212/729
1    |  9017/3168     −355/33    46732/5247     49/176   −5103/18656
1    |    35/384         0         500/1113    125/192   −2187/6784     11/84	
--------------------------------------------------------------------------------------
     |  5179/57600       0        7571/16695   393/640  −92097/339200  187/2100  1/40
     |    35/  384       0         500/1113    125/192   −2187/6784     11/84    0

0  |
c2 | c21
c3 | c31 c32
c4 | c41 c42 c43
c5 | c51 c52 c53 c54
c6 | c61 c62 c63 c64 c65 
c7 | c71 c72 c73 c74 c75 c76
----------------------
   | c81 c82 c83 c84 c85 c86 c87
   | d81 d82 d83 d84 d85 d86 d87

Функции h2-h8:

h2(mp, dt, c21, p1') {
    return mp + dt * (c21 * p1')
}
h3(mp, dt, c31, p1', c32, p2') {
    return mp + dt * (c31 * p1' + c32 * p2')
}
h4(mp, dt, c41, p1', c42, p2', c43, p3') {
    return mp + dt * (c41 * p1' + c42 * p2' + c43 * p3')
}
h5(mp, dt, c51, p1', c52, p2', c53, p3', c54, p4') {
    return mp + dt * (c51 * p1' + c52 * p2' + c53 * p3' + c54 * p4')
}
h6(mp, dt, c61, p1', c62, p2', c63, p3', c64, p4', c65, p5') {
    return mp + dt * (c61 * p1' + c62 * p2' + c63 * p3' + c64 * p4' + c65 * p5')
}
h7(mp, dt, c71, p1', c72, p2', c73, p3', c74, p4', c75, p5', c76, p6') {
    return mp + dt * (c71 * p1' + c72 * p2' + c73 * p3' + c74 * p4' + c75 * p5' + c76 * p6')
}
h8(mp, dt, c81, p1', c82, p2', c83, p3', c84, p4', c85, p5', c86, p6', c87, p7') {
    return mp + dt * (c81 * p1' + c82 * p2' + c83 * p3' + c84 * p4' + c85 * p5' + c86 * p6' + c87 * p7')
}

Метод ДП:

dormandPrince(mp, mt, dt) {
    p1' = g(   mp                                                                 , mt          )
    p2' = g(h2(mp, dt, c21, p1'                                                  ), mt + c2 * dt)
    p3' = g(h3(mp, dt, c31, p1', c32, p2'                                        ), mt + c3 * dt)
    p4' = g(h4(mp, dt, c41, p1', c42, p2', c43, p3'                              ), mt + c4 * dt)
    p5' = g(h5(mp, dt, c51, p1', c52, p2', c53, p3', c54, p4'                    ), mt + c5 * dt)
    p6' = g(h6(mp, dt, c61, p1', c62, p2', c63, p3', c64, p4', c65, p5'          ), mt + c6 * dt)
    p7' = g(h7(mp, dt, c71, p1', c72, p2', c73, p3', c74, p4', c75, p5', c76, p6'), mt + c7 * dt)

    //считаем два решения
    pp1 = h8(mp, dt, c81, p1', c82, p2', c83, p3', c84, p4', c85, p5', c86, p6', c87, p7')
    pp2 = h8(mp, dt, d81, p1', d82, p2', d83, p3', d84, p4', d85, p5', d86, p6', d87, p7')

    return (pp1, pp2);
}

Продолжение следует...

Верхний уровень алгоритма

ДП возвращает два вектора pp1 - решение, pp2 нужен только для расчёта ошибки:

stepAndError(p0, t, ∆t) {
    (pp1, pp2) = dormandPrince(p0, t, ∆t) 
    return (pp1, |pp1 - pp2|)
}

Всё.

Вот это уже разговор. :-) Спасибо, Станислав. Какие у меня мысли перед кодированием:

  • pp1 и pp2 записываем в старые p1 и p2,
  • константы cjk хорошие. Мы будем хранить их не в регистрах, а в коде. С этого, пожалуй, и начну. В регистры же можно записать адреса подпрограмм, их вычисляющих.
  • Следующим этапом идёт кодирование hn(), объединив их все в одну подпрограмму — это самое сложное.
    • В R0 она будет получать, сколько производных ей складывать.
    • Результат она будет готовить сразу в q на вход g().
    • Вполне возможно, что каждую из компонент каждой производной мы сможем посчитать с точностью 14 разрядов и обрезание произойдёт лишь при записи результата в регистр.
  • g() пока возьмём старую. Хотя она могла бы результат записывать с помощью косвенной адресации… тут подумаю. Похоже, что придётся опять перетасовывать производные и облагораживать Рунге-Кутта.
  • девять вызов из основной подпрограммы метода писать будет просто, при правильном подходе там векторов практически не будет видно.

> Следующим этапом идёт кодирование hn(), объединив их все в одну подпрограмму — это самое сложное.
> Результат она будет готовить сразу в q на вход g().

Вот возможный дизайн:

q = 0 ; вектор, куда будут складываться производные
x = адрес конца вектора [c71, c72, c73, c74, c75, c76]
ПП h7

h8:
   q = q + [x--] * p7'
h7:
   q = q + [x--] * p6' ; извлекаем регистр по адресу x, уменьшаем x на один
h6:
   q = q + [x--] * p5'
h5:
   q = q + [x--] * p4'
h4:
   q = q + [x--] * p3'
h3:
   q = q + [x--] * p2'
h2:
   q = q + [x--] * p1'
   q = mp + dt * q
   В/O

Операцию q = q + c * p' можно закодировать подрограммой

; заталкиваем c
; заталкиваем адрес p'
ПП inc

:inc
    q.x = q.x + c * p'.x ; адреса можно перебрать автоинкрементом
...
    В/O

Получится довольно компактно.

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

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

— не люблю ими жертвовать. Пока расписал без циклов, оказалось не так сложно.

Когда вы дадите добро на математику, начну изменять номера регистров и адресация станет простой. Сделаем циклы. Все производные можно вообще выселить за R100 — насколько я понимаю, вне h() к ним никто никогда обращаться не будет. Ну ещё новая g() туда результат запишет — опять же, с помощью косвенной адресации.

Нужны тесты, чтобы узнать чем вы пожертвовали и что выиграли.

Наврали здесь:

c2: 0 RTN ; 0
c3: 5 F1/x RTN ; 1/5

Исправил. Загрузил версию 0.14

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

В результате оптимизации мы можем прийти к чему-то похожему. Возможно, действительно придётся написать макрогенератор — или воспользоваться чем-то вроде m4. Возможности .EQU не позволяют вставлять коэффициенты прямо в код, как inline. Но всё это позже.

Пусть процедура stepAndError написана. В ней нет условной логики. При любых входных данных она совершает одну и туже последовательность действий. Для этих действий нужны шесть аргументов (константы я игнорирую), результат - шесть значений. Раскроем эту последовательность (например автоматически) в простыню из команд. Можно написать программу, которая автоматически оптимизирует (минимизирует) перемещение данных между регистрами и использование самих регистров. То что получится, будет большой кашей из инструкций, но кашей работающей правильно и быстро (возможно, самой быстрой из возможных). Видимо, объём программы не важен и такой подход имеет право на существование.
Если вам это интересно, то я могу попробовать решить задачу автоматической оптимизации.

Я реализовал макросы с помощью косвенной адресации. :-) Где нужно, например, умножить на 5179/57600 /* коэффициент c81 */, я пишу КПП7.

При этом, разумеется, пришлось записать эту подпрограмму отдельно, закончив В/О и заранее занеся её адрес в R7.

Сейчас 14-разрядные расчёты каждой из производных происходят с 14-разрядными коэффициентами, но 12-разрядными значениями предыдущих производных, что немного снижает точность. После отладки попробую привести в порядок номера регистров и свернуть каждый из h2…h8 в 4-х проходовый цикл. Тогда можно будет перенести коэффициенты внутрь программы и сделать 14-разрядную точность не только для расчёта производных, но также для их хранения.

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

Ловите версию 0.15 с реализацией Дорманд-Принса. Я её ещё не запускал.
Когда отладим, будем добавлять дополнительные блоки и оптимизировать.
http://the-hacker.ru/2011/lun-103.mkl
http://the-hacker.ru/2011/lun-103.mkp
http://the-hacker.ru/2011/lun-103.htm

P.S. Программа занимает 24 страницы, почти четверть всей программной памяти МК-161. Мне не пришлось использовать регистры R100 и старше, но регистры до R99 включительно практически исчерпаны. Помимо трёх неиспользованных регистров (R15, R68 и R69) удалось сохранить младший регистровый НЗ на циклы и автодекремент/автоинкремент (R0…R6) — который теперь будем свободно тратить на оптимизацию.

Отлично летает! Вот тестовый полёт на "сетевом" "лунолёте".

В программе я убрал скорости на старте, залил полный бак (3500), установил метод ДП (2), ошибки уменьшил до 10^-9.

Команды:

r/s
84 r/s 2100 r/s 74 r/s
0 r/s 0 r/s 2035 r/s
90 r/s 46 r/s 2 r/s
0 r/s 0 r/s 7266,1717 r/s
0 r/s 0 r/s 0,0000625 r/s

Последняя команда вызвана тем, что эмулятор не принимает больше 8 цифр в одной константе?!
(А можно чтобы точка в числах читалась? А почему она тогда печатается?)

Результаты:

vy -0.01634821061063
vx  1619.1402936749
y   133511.2528215
x   3195583.1328723

х я считал так: PRM26 FPI 2 * - PRM42 *

Совпадение отличное!

Просьба к владельцам калькуляторов: сравните производительность методов, пожалуйста.

Ура! А ведь по этому поводу можно и сходить за черноголовской бутылочкой «Тархуна» или «Байкала»! :-)) Записал ваши начальные данные в программу, установил Дорманд-Принс методом по умолчанию, поднял точность до 1 ВП –07, откомпилировал и успешно повторил ваш первый полёт. Новый «Лунолёт» версии 0.17 с этими настройками и другими незначительными изменениями уже лежит на сайте.

Чтож. Теперь начинаю выносить регистры производных за барьер R100 — давая им такие номера, чтобы упростить hn() и подготовить их к переходу на 14 разрядов. Огромное вам спасибо за уже проделанную работу и надеюсь, что вы меня не бросите до релиза. Понадобится и тестирование оптимизированных версий, и математика для посадки. Далее, ответы на ваши вопросы по эмулятору.

1. Экран МК-161 отображает только 8 цифр мантиссы RX, зато крупным шрифтом. Соответственно ввести с клавиатуры можно только эти 8 цифр и в эмуляторе есть специальный код, блокирующий остальные цифры для совместимости с ЭКВМ.

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

7266,1717 В↑ 0,0000625 +

2. Во входном языке MK.EXE используется исключительно десятичная запятая, точка зарезервирована для начала псевдокоманд. ВК-6 также строго следует этой дисциплине, чтобы не приучать владельца ЭКВМ к ошибкам. Возможно, я сделаю точку первым символом для команд отладчика, чтобы этот вопрос не возникал. (записал в todo.txt)

3. Вывод стека происходит сейчас с помощью стандартной функции printf() и был написан, как нечто временное — пока кто-нибудь не сделает настоящий вывод в графическом интерфейсе. Поскольку графический интерфейс магически не возник и временное стало постоянным, записал ваше предложение в todo.txt

> Вообще-то на экране МК-161 есть место и для девятой цифры, но даже если мы все дружно попросим Михаил Борисовича её задействовать, погоды она не сделает.

Давайте попросим Михаила Борисовча сделать ввод 12 цифр, их скролирование на экране и возможность забоя введённых цифр (и ВП и порядок и +/- включить в туже схему с возможостью забоя). Забой не нужен в виде команды программы, но удобный буфер ввода с минимальным редактированием (забоем) очень желателен.

> 2. Во входном языке MK.EXE используется исключительно десятичная запятая, точка зарезервирована для начала псевдокоманд. ВК-6 также строго следует этой дисциплине, чтобы не приучать владельца ЭКВМ к ошибкам. Возможно, я сделаю точку первым символом для команд отладчика, чтобы этот вопрос не возникал. (записал в todo.txt)

Когда-то точка и запятая расколола комитет по стандартизации АЛГОЛА-68.
Обработка и точки и запятой внутри числового литерала, на мой взгляд, не приводит к противоречиям во входном языке эмулятора или компилятора. Но если такова политика, то я сожалею.

> Поскольку графический интерфейс магически не возник и временное стало постоянным, записал ваше предложение в todo.txt

Не воспринимайте текстовый интерфейс как временное решение. Он хорошо подходит для организации взаимодействия эмулятора с другими программами на компьютере. Микропример: для тестов я запускал эмулятор и в его консоль копировал последовательность команд (в одну операцию). Это серъёзно ускорило тестирование. Графический интерфейс заставил бы меня вводить команды по одной, что замедлило бы тестирование и сделало бы его менее качественным (баг с мусором помните?). А текстовый интерфейс позволит при случае построить автоматические тесты, когда компьютер готовит данные вводит их в эмулятор, считывает ответы, проверяет их правильность.

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

Кстати, возможен ли сейчас такой тестовый сценарий с участием калькулятора и компьютера?


> Давайте попросим Михаила Борисовча сделать ввод 12 цифр, их скролирование на экране
> и возможность забоя введённых цифр (и ВП и порядок и +/- включить в туже схему
> с возможостью забоя).

Попробуйте. Это серьёзная переделка режима F АВТ, но я — только за. Тем более, что клавиша ← на клавиатуре ЭКВМ есть. Забой мне несколько раз требовался при работе на МК-161, а ввод длинных цифровых «хвостов» постоянная головная боль. Я ввёл это в ВК-6 с большим сожалением и лишь для того, чтобы отлавливать эти ошибки в программах для ЭКВМ. Ограничить точность регистров, как видите, только сейчас решился.

> Обработка и точки и запятой внутри числового литерала,

Сейчас в MK.EXE нет числовых литералов. 12 это две команды с кодом 01 и 02. Точно также числа воспринимает и ВК-6.

Введение литералов и изменение паржинга псевдооператоров просто, чтобы ввести в отечественный ЯМК американскую точку, как единственную дублирующуюся мнемонику ЯМК — не лучше ли это время потратить на добавление чего-нибудь действительно полезного? Тем более, что даже американцы сейчас вводят запятую в локализацию своих операционных систем.

> Не воспринимайте текстовый интерфейс как временное решение.

Да, я это уже понял. :(

> Кстати, возможен ли сейчас такой тестовый сценарий с участием калькулятора и компьютера?

Когда я работал под Windows, я заливал свежий откомпилированный код в МК-161 и тут же его проверял. Вот и вся автоматизация того времени.

Входной поток ВК-6, разумеется, можно подавать из файла или другой программы.

>> Давайте попросим Михаила Борисовча сделать ввод 12 цифр, их скролирование на экране
>> и возможность забоя введённых цифр (и ВП и порядок и +/- включить в туже схему
>> с возможостью забоя).

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

Вообще говоря, моё личное мнение здесь роли не играет.

Основной аргумент против ввода 12 и более символов в RX при сохранении вывода 8 цифр - это усложнение набора (с использованием сдвига строки или переключения размера шрифта), чреватое возникновением ошибок при вводе. Вероятно, следует напомнить, что режим индикации с выводом 14 цифр в регистрах Y, Z и T изначально был недокументирован и предназначался только для отладки.

Против вывода более 8 символов в RX - недостаток места для вывода шрифтом 1; предназначение ЭКВМ для инженерных, в первую очередь, расчётов; совместимость с ПМК.

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

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

> Основной аргумент против ввода 12 и более символов в RX при сохранении вывода 8 цифр - это усложнение набора (с использованием сдвига строки или переключения размера шрифта), чреватое возникновением ошибок при вводе.
Добавление забоя сохраняет обратную совместимость с точки зрения пользователя: любая последовательность допустимых действий сейчас сохранит свой смысл после ввода забоя.

> Вероятно, следует напомнить, что режим индикации с выводом 14 цифр в регистрах Y, Z и T изначально был недокументирован и предназначался только для отладки.
Сейчас это важный элемент дизайна, который облегчает жизнь пользователю - вместо одного числа видна картина целиком. Это важно и при вводе и при выводе.

> Против вывода более 8 символов в RX - недостаток места для вывода шрифтом 1; предназначение ЭКВМ для инженерных, в первую очередь, расчётов; совместимость с ПМК.
Недостаток места: уменьшайте шрифт после ввода 9 цифры, или скролируйте.
Инженерные расчёты: Разрядов ещё никому не было достаточно.
Совместимость с ПМК: в данном случае это вредная совместимость, о потере которой никто не пожалеет.

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


> Основной аргумент против ввода 12 и более символов в RX при сохранении
> вывода 8 цифр - это усложнение набора (с использованием сдвига строки
> или переключения размера шрифта), чреватое возникновением ошибок при вводе.

Тут как раз будет упрощение набора. Люди привыкли к изменению размера шрифта (и появлению второй строчки) при вводе номера телефона в мобильных аппаратах.

А вот разбитие 14-разрядного числа при вводе на пару чисел и просчитывание порядка, чтобы «хвост» лёг, куда нужно — тут ошибки действительно возникают. Не говоря уже о том, что как-то унизительно пользователя заставлять это проделывать, в XXI веке.

> Вероятно, следует напомнить, что режим индикации с выводом 14 цифр
> в регистрах Y, Z и T изначально был недокументирован и предназначался только для отладки.

Режим с выводом 8 знаков мелким шрифтом, с огромным пробелом при 14 вычисленных знаках — это варварство! Несколько лет назад попросил вывод 14-разрядных (или хотя бы 12-разрядных) чисел в RY, RZ и RT в естественном виде. Чтобы 100500 отображалось, как 100500, а не 1,005 05. Но до сих пор приходится «запятую» переносить самостоятельно, в уме.

> Ввод команд в режиме АВТ обрабатывается тем же способом, что и в режиме ПРГ.

В режиме ПРГ, кстати, можно нажать клавишу «назад» и исправить цифру вводимого числа. А вот в режиме АВТ такая возможность отключена. Можно даже, чтобы «забой» работал лишь на мантиссу, порядок ещё в ПМК можно было перезабить. Возможность исправление ошибок при вводе существенно поможет их отсутствию — то есть декларированной цели.

Про дурную совместимость соглашусь со Станиславом. Я не буду жалеть, если десять команд "1" введут 10-значное число, а не 8-значное. Всё равно никто в здравом уме таких «программ» для ПМК не писал. А вот для ЭКВМ мне приходится проверять, укладываются ли коэффициенты Дорманд-Принс в прокрустово ложе из 8 цифр или нет — оказалось, там лишь 6 знаков и я вздохнул с облегчением. Менее же опытный программист однажды нарвётся здесь на ошибку.

Отвечу сразу на оба сообщения.

Доводы "за" мне тоже известны. Кроме того, моё личное мнение может не всегда совпадать с мнением конторы. ;-)

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

Коэффициенты d81 d82 d83 d84 d85 d86 совпадают с c71 c72 c73 c74 c75 c76, а d87 вообще равен 0.

Означает ли это, что pp2 равно первому аргументу g() в вычислении p7' и я впустую гоняю МК-161, рассчитывая одно и тоже второй раз?

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

Выложил версию 0.16 с этой оптимизацией. Программа похудела на 1 страницу.

Жду-недождусь, когда освободимся от балласта Эйлера и Рунге-Кутта. :-)

Программа летает. Результаты похожи на "сетевые". Разница объяснима, так как точность в программе 10^-7, а в сети 10^-9.

Но есть и проблема.
Я загрузил программу и выполнил команды:

r/s
84 r/s 2100 r/s 74 r/s
0 r/s 0 r/s 2035 r/s
90 r/s 46 r/s 2 r/s
0 r/s 0 r/s 7266,1717 r/s
0 r/s 0 r/s 0,0000625 r/s
PRM26 FPI 2 * - PRM42 *

Эта последовательность считает координату X после витка на орбите.
Результат: X 3195583.1311796

Потом я строку

                0,000001 PM43           ; dt_min

заменил на

                6 +/- F10^x PM43        ; dt_min

и проделал тот же тест.
Результат: X 3195583.1319331

Отличаются последние четыре знака. В чём дело?

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

http://narod.ru/disk/30542552001/lun-103.zip.html

Баг исправлен

В версии 0.19 я существенно перетряхнул регистры и сжал все функции h() в одну hn(). Программа теперь занимает всего 20 страниц, пятую часть памяти.

Вместе с тем все вычисления внутри всех трёх методов теперь проводятся с полной точностью стека (14 разрядов). Для хранения компонент производных и других промежуточных вычислений используется по два регистра. Когда код будет отлажен, можно будет вынести 14 разрядов в основной цикл и легко увеличить точность на два порядка (например, 1 ВП –09 вместо 1 ВП –07). Тогда скорости на экране будут содержать все 14 разрядов и плюс два разряда к высоте. :-)

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

Качать здесь:

http://narod.ru/disk/30549540001/lun-103.zip.html

У меня на эмуляторе программа выполнила первый манёвр (84 r/s 2100 r/s 74 r/s) правильно, но очень долго (∆t = t_min ?) и зависла на втором (0 r/s 0 r/s 2035 r/s). :(

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

> Новая версия часто ошибается с прогнозом шага.

Это лишь на ДП или на РК и Э тоже? В математике я не менял ничего и прогнозы должны незначительно отличаться от того, что выдавала предыдущая версия.

Для РК поведение такое же.

Вы выбрали неверное направление. 14 знаков против 12 не дадут большого выигрыша в точности. Прежде чем, оптимизировать точность я бы провёл эксперименты на ВК-6, зарезая точность всех вычислений до 12 и до 14 знаков. Такой тест дал бы основания для работ по улучшению точности.
Сейчас важнее производительность - коэффициенты в регистрах, быстрое вычисление линейных комбинаций, минимум лишних перемещений между регистрами. Если оценки производительности ДП будут неутешительными, его вообще придётся выкинуть (какая уж тут точность!).

Версию 0.18 я у себя сохранил — тестировать и развивать 12-разрядную версию можно параллельно.

Просто числодробилка МК-161 может выдать максимум 14 разрядов и если быстродействие позволит (объём памяти позволяет), хотелось бы использовать данный ресурс полностью. Из обоих версий можно выкинуть как ДП, так и РК.

14-разрядная версия 0.19 выложена на хостинге:
http://the-hacker.ru/lun-103.html

> эксперименты на ВК-6, зарезая точность всех вычислений до 12 и до 14 знаков.

Если под «всеми вычислениями» вы имеете в виду запись чисел в регистры, то инструкции по модификации ВК-6 я выложил. Хоть 12, хоть 14 — сколько угодно знаков.

В версии 0.19 точность вычислений в стеке не изменилась. Просто я записываю 14-разрядный ответ не в один 12-разрядный регистр, а в пару по этому методу. Само по себе это точности решений не даст, но после внесения всех изменений и отладки позволит находить два дополнительных знака решений, то есть уменьшать запрашиваемую погрешность (R34 и R35) на два порядка по сравнению с 12-разрядной версией.

В новой версии баг остался и на РК и на ДП.

Большая таблица ДП может привести к медленному счёту, но есть методы четвёртого порядка с меньшими таблицами: Cash-Karp и Runge-Kutta-Fehlberg.

Это та же самая версия, что и на народе — с нелокализованной ошибкой. А Эйлер работает?

Вот эта строка вызывает подозрения:

		PRM36 PRM37 * PRM50 PRM51 * - Fx^2		; (p1.x * p1.y - p2.x * p2.y)^2

При РК регистры 36 и 50 могут содержать значение 257631 (радианы). Значения в регистрах равны.
При ДП регистр 36 и 50 могут содержать 658.08126065432 (радианы). Значения в регистрах равны.
Похоже, что какой-то фрагмент затирает значения в регистрах.
Эйлер ошибок, кажется, не делает.

Правильно ли я понимаю, что все вызовы KGSBe выдают чепуху в q.x = R90+R91 — откуда она потом разносится в p1.x и p2.x (указанные регистры)? В этом случае ошибка РК лезет из этого фрагмента LinEx (1072):

                79 M4                                           ; R80,81 = p'x'
                KRM4 KRM4 + PRM67 * PRM20 +                     ; q.x  := mp.x  + mdt * p’.x’
                PM90 PRM90 - PM91                               ; R90,91 = q.x

Ещё можно проверить, когда Fizika начинает врать при расчёте координаты x (1568):

                PRM94 PRM95 + PRM92 PRM93 + /           ; pn'.x' := q.vx / q.y
                KM8 KRM8 - KM4                          ; записать в двойной регистр pn'.x'
 

Параллельно с этой «трудной» 14-разрядной версией (я предупреждал, что перетасовка регистров будет наиболее сложным моментом) я выложил новую 12-разрядную версию.

В ней я учёл ваши советы по повышению быстродействия — занёс коэффициенты в регистры, вставил h() как inline-код в метод, оптимизировал случаи когда cxy=0 и т.д. Также я создал новую тему и предлагаю работать над программой в ней. Не знаю, как у вас, а у меня Друпал глючит, когда комментариев так много, что они переходят на другую страницу и искать ваши ответы здесь превращается в отдельную компьютерную игру.

Сравнил шаги, которые выбирают РК и ДП при экстраполяции. Они почти одинаковы. Следовательно точность обоих методов тоже примерно одинакова (что, в частности, означает, что оба метода реализованы одинаково правильно).

Замечу, что я недавно уповал на то, что ДП даст больший шаг. Это было глупостью с моей стороны. Методы одного порядка должны давать одинаковые оценки ошибки.

С нетерпением жду сравнения производительности методов на калькуляторе.

P.S. Было бы удобно, если эмулятор умел оценивать время работы программы на калькуляторе.

Тогда в чём выгода более сложного ДП, если методы одного порядка дают одинаковые точность, ошибки и шаг? Ещё не поздно его убрать.

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

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

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

Если у нас шаги одинаковые с РК, значит что-то не так с реализацией ДП.

> Я понял теорию так, что порядок метода даёт коэффициент при вычислении шага.

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

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

Таблица обеспечивает порядок метода и, что важнее, его корректность. Я пробовал менять числители коэффициентов - лунолёт бодро прилетал не туда вообще. Разработка таблиц дело трудное - все таблицы именные.

> Если у нас шаги одинаковые с РК, значит что-то не так с реализацией ДП.

Всё с ДП так. РК и ДП дают почти одинаковую оценки ошибки. Зато ДП считает более точно (для продолжения решения берётся ветка пятого порядка). ДП вычисляет функцию g() 8 раз (7 с учётом вашего замечания), а РК - 12 раз. ДП проигрывает РК в количестве линейных комбинацй.

ДП начнёт выигрывать в будущем, когда в уравнения попадут вместе и Луна и Земля с их вращениями ("Лунолёт-4").
Чем сложнее g(), тем выгоднее ДП.

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

Спасибо за объяснения.

Любопытно, рассматривали ли вы увеличение разрядности мантиссы первых 15 регистров до 14 десятичных разрядов. Это на два порядка поднимет точность простеньких программ — навроде тех, что в справочнике проф. Дьяконова.

> "навроде тех, что в справочнике проф. Дьяконова."

Типа расчёта числа витков катушки индуктивности? До 12-го или 14-го разряда. ;-)

> "увеличение разрядности мантиссы первых 15 регистров до 14 десятичных разрядов."

Обращение ко всем регистрам происходит одинаковым образом. Различие по группам 0-15, 0-99, 0-9999 существует только при обработке кода команд, при формировании адреса регистра.

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

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

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

Но вот вы делаете сейчас программу расчёта траектории КА, с вычислениями до 12-14 знака. Не буду даже упоминать, что отработать манёвр с такой точностью реальный двигатель не сможет, масса аппарата и прочие параметры известны лишь приблизительно, небесные тела несферические и имеют неоднородное поле и даже гравитационная постоянная определена земной наукой только с 3-4 знаками:

"В единицах СИ рекомендованное Комитетом данных для науки и техники (CODATA)[3] на 2008 год значение было G = 6,67428(67)·10−11 м³·с−2·кг−1, или Н·м²·кг−2. в 2010 году значение было исправлено на: G = 6,67384(80)·10−11 м³·с−2·кг−1, или Н·м²·кг−2. В октябре 2010 в журнале Physical Review Letters появилась статья[4], предлагающая уточнённое значение 6,67234(14), ... Пересмотр величины G, произошедший в период с 1986 г. по 2008 г., был вызван исследованиями неупругости нитей подвесок в крутильных весах[5]." (википедия)

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

P.S. С праздником!

Вот примерно поэтому я и не использую в программе ни g, ни G — перейдя на более точно известные GE и μ. :-) Но даже независимо от значения постоянных, модель должна подчиняться законам классической механики — а не разрядной сетки ЭКВМ. Масса корабля может быть недостаточно хорошо известна пилоту, но она идеально точно известна как притягивающей её планете, так и работающему двигателю. Наша задача — моделировать эти физические явления, а не хрупкую психику пилота.

Если бы мы ограничили модель природы человеческим невежеством и задали точностью вычислений на секунду ваши три знака, то за 17 минут такая модель начала бы ошибаться в разы. Чтобы летать всего сутки, запас точности должен составлять пять знаков. Миссия же к Марсу занимает годы, это ещё 3-4 знака. Раскладки же про хранение декартовых координат с таким запасом точности, учитывающие размер Солнечной системы, на этом форуме приводились.

Например, мне хотелось бы не просто, скажем, знать положение Сатурна относительно Солнца с точностью до миллиметра. Мне хотелось бы промоделировать хотя бы один его оборот вокруг Солнца в простых и надёжных декартовых координатах. Сколько дополнительных знаков нужно к разрядной сетке, вам посчитать или справитесь? Или я на слишком большие задачи замахиваюсь, такое считать людям на калькуляторах не позволено? А я ведь ещё о полётах к звёздам задумываюсь, и о других галактиках… :-)

Изменяя подпрограмму Fizika можно учитывать столько факторов, сколько сможет придумать фантазия. Но хочется, чтобы в пределах каждой из задач владелец ЭКВМ получал бы точное и чистое представление о том, как действует каждый из учтённых факторов — а не воздействие шумов от ограниченной разрядности.

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

Вас тоже с наступившем.

> На мой взгляд, вы просто пытаетесь перевести тему с обсуждения ЭКВМ на незавершённый проект, пользуясь его открытостью

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

> Поэтому без вашего участия не выяснить, что мешает сделать при обращении к десятичным регистрам ветвление на группу 0-15, сравняв их разрядность с разрядностью стека.

Технически, как понимаете, ничего не мешает. Сотня байт в памяти нашлась бы, скорее всего.

Если для марсианской миссии потребуются ЭКВМ - так и быть, сделаем. ;-)

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

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

> Нет. Просто пытаюсь в очередной раз обратить Ваше внимание на нецелевое использование ЭКВМ - она пригодна скорее для управления реальным Лунолётом, чем для подобных вычислений.

Хорошая математическая библиотека могла бы раширить (потенциальную) аудиторию пользователей. Методы интегрирования систем дифференциальных уравнений - часть такой библиотеки. Качественные методы - конкуретное преимущество.

Понятно, что такие вещи как расширение рынка звучат очень туманно. И в любом случает понадобится масса усилий и воли, чтобы хотя бы попытаться раширить аудиторию. И наверняка попытка провалится. С другой стороны почти все начинания, которые сейчас успешны, начинались с множества провалившихся попыток. Поразительно, что вам вообще удалось сделать калькулятор. Надо продолжать делать глупости :) с его применением. Одна из них может стать решающей и идея выстрелит (или проживёт на несколько лет дольше).

Компьютеры (калькуляторы) как устройства для обработки данных (вычислений) - до сих пор ниша, которая не заполнена. ЛЬвиная доля цифровой техники играет роль абонентских усройств - устройств передачи данных.

Про управляющие устройства я и не говорю - это, наверно, ваш рынок и есть.

Устройства для обработки данных существуют - это компьютеры и калькуляторы. Не все ещё персональные компьютеры дошли до стадии абонентских устройств, хотя "невидимая рука" их туда стремительно впихивает.

Управляющие машины тоже существуют во множестве. Но большая часть промконтроллеров стоит довольно дорого и имеет специфические особенности.

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

Я надеюсь, что применение МК-161 для расчёта космических путешествий станет столь же важным и начнёт оказывать влияние на сбыт и разработку ЭКВМ.

Здесь не только вечная мечта о космосе. Мой непростой жизненный опыт показывает, что ум хэкерам имеет смысл развивать именно на космических приложениях. Альтернативный криминальный земной путь завязан на такое количество конфликтующих коммерческих и милитаристских интересов, что юный хэкер просто не доживёт до какого-либо существенного уровня — по крайней мере сохранив свою творческую свободу, независимость и интерес к искусству.

Когда ты вычисляешь тонкости движения Луны, корабля в криволинейных координатах или обрабатываешь снятый на веб-камеру видеоролик Юпитера, это вызывает меньше страстей и агрессивной животной реакции «нас посчитали», чем если ты изучаешь систему защиты, реализованную корпорацией Сони в PlayStation 3.

Для дальнейшей мотивации дискуссии :)
http://ru.euronews.net/newswires/1207731-newswire/

"Сотрудники Системы дальней космической связи НАСА послали команду на переключение 4 ноября. Но только 7 ноября они получили ответ “Вояджера”, подтверждающий, что он воспринял команду и переключился на резервную систему контроля ориентации аппарата. Столько дней и часов потребовалось радиосигналу, мчащемуся со скоростью света /300 тыс км в сек/, чтобы догнать аппарат и принести на Землю его ответ. Выполненное переключение даст возможность экономить примерно 12 ватт энергии, что, по словам представителей НАСА, позволит “Вояджеру” прожить еще одно десятилетие."

---------------------------
Истина где-то рядом
www.litres.ru/vitaliy-samurov/dozvonitsya-do-devy/

Сугубо для сравнения отправил Вам в ЖЖ сообщение с блоком подпрограмм, которые с некоторым приближением вычисляют орбиту Луны. Примерно так же выглядят и расчёты траектории КА, в дополнение к обычной классике по ТАУ. ;-)

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

Но сейчас, при моделировании межпланетных полётов, нам придётся идти на множество хитростей — даже если МК-161 потянет по быстродействию все 14 разрядов. Просто не надо делать вид, что задач, требующих высокой разрядности, не существует. Существуют такие задачи, особенно в численном моделировании. Но пока по одёжке протягиваем ножки. И на счётах считали. Но, как уже было здесь сказано, разрядности много не бывает.

Задачи такие, конечно же, существуют.

Более того, любая носогрейка из "века четырёхгигабайтных флешек", в принципе способна их решать. Только вот реальной потребности в таких расчётах нет, к сожалению.

Михаил Борисович, спасибо вам и за 12 разрядов. :-) Уже стало заметно легче, чем на МК-61.

Я ведь просто поинтересовался про запас точности традиционных регистров. Не настаиваю и понимаю, что всем не угодишь. Для меня это не первостепенная важность. В моей задаче 15 регистров погоды всё равно не сделают, а хранить 14 разрядов в двух регистрах я и сам научился.

Спасибо не мне, но передам. ;-)

Понятно, что совершенствоваться в каждую из сторон можно до бесконечности. Число разрядов в регистрах ЭКВМ - такой же компромисс, как и многое другое.

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

Время же и пространство бесконечны и неисчерпаемы, ограничены лишь наши вычислительные возможности. :-)

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

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

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

Так куда уж проще использовать Питон на любой мало-мальски современной железяке? :)

Python 2.7.1 (r271:86832, Nov 27 2010, 18:30:46) [MSC v.1500 32 bit (Intel)] on win32
Type "copyright", "credits" or "license()" for more information.
>>> earth=(360/365.0)*60
>>> mars=earth/1.881
>>> diff=earth-mars
>>> print diff
27.7171134561
>>>

---------------------------
Истина где-то рядом
www.litres.ru/vitaliy-samurov/dozvonitsya-do-devy/

Это тоже имеет право на существование, разумеется.

Относительно же простоты - вопрос спорный. У Питона есть масса недостатков, его главное преимущество - открытость.

Но давайте эту тему холиварами не засорять - она и так большая. ;-)

> Пока сделал Рунге-Кутта методом по умолчанию и поставил 10–6, половину 12-разрядной сетки, а там посмотрим. У меня даже в этом значении сомнения, т.к. радиус Луны много отъедает точности от координаты y, а ведь по ней рассчитывается гравитация и центробежка.

В данном случае, на точность расчётов влияет не абсолютная погрешность (10^-5 в случае поверхности Луны), а относительная (10^-11). Специфика чисел с плавающей запятой такова, что относительная погрешность не зависит от порядка числа - она примерно одинакова почти во всём диапазоне значений. Расчёт ускорения при этом страдает мало.

Вы правы, про ускорения я несколько погорячился. Жизнь не так плоха. Возможно, 12-разрядной сетки МК-161 хватит на координату y даже без трюков. Посмотрим. Хорошо, что у нас всегда есть возможность сравнить с компьютерными вычислениями. Конечно, если вы поправите на сайте гравитационную константу и радиус Луны под новые данные.

Вот немножко теории на этот счёт:
http://lunolets.livejournal.com/1012.html

Вот комплект данных, которые надо ввести (вставить) в окно "сетевого" "лунолёта" для совершения тестового манёвра:

R 1737400
g 1.624216906405779
v 100 u 50
h 0 x 0
m 100
--2
45 90 5

Совпадение результатов с вашим последним вариантом я проверил

Для получения из ВК-6 точного эмулятора достаточно написать функцию real Shrink12(real), округляющую числа до 12 десятичных разрядов мантиссы и добавить её в файл reg.c

Самое простое, вероятно, это сделать через sprintf/sscanf (я этот код не отлаживал, просто прикинул):

long double Shrink12(long double x)
{ static char s[40];  static long double y;
  sprintf(s, "%.12LG", x);
  sscanf(s, "%LG", &y);
  return y;
}

Потом в этом же файле reg.c заменить ldR[r] = RX; на ldR[r] = Shrink12(RX); в следующих трёх функциях:
static void StoreR(int r)
void M(int r)
void PM(int bcd)

Ну и откомпилировать его, конечно. Проверить усечение разрядов можно, например так: FPI M9 RM9 —

Где мне почитать про форматы хранения чисел в стеке и регистрах калькулятора?

http://mk.semico.ru/mk161re.htm#p2
12 десятичных знаков мантиссы + 2 порядка (в регистрах)
14 десятичных знаков мантиссы + 2 порядка (в стеке)

Будете ли вы встраивать усечение в основную версию ВК-6?
Будете ли усекать регистровые операции?

Я написал себе в todo, чтобы сделать такую опцию в командной строке. Но код, который я написал выше, работает.

У меня нет эмулятора с собой. На работе я работаю :).
Я попробовал на сайте ввести команду 1 90 5 и получил в точности вам результат для РК.

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

dt_min — наша гарантия от зацикливания. Каждый шаг будет продвигать нас к концу манёвра хотя бы на чуть-чуть. Да, при выходе на dt_min точность уменьшится. Зато машинка не зациклится и даже если считает долго, результата дождаться можно.

min_dt - форма обмана. Я предлагаю останавливать программу с ошибкой если dt обнулился.

1. dt никогда не обнулится если заданы вменяемые параметры (скорости, ускорения и расстояние от центра тяготения)
2. Если min_dt сработает хоть раз, гарантии на точность заканчиваются.

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

Это не единственная форма обмана. Например, у нас ещё разрядная сетка ограничена.

Опасно не только зануление ∆t. Существуют бесконечные ряды, сумма которых — конечное число. Чтобы доказать, что наша реализация их не выдаст, надо затратить довольно большие усилия.

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

Да есть ещё округления. И с ними трудно бороться - хотя можно оценить их вклад в ошибки. Я про то что у нас есть обман, который легко убрать.

Наличие большого зла в мире не есть повод не бороться со злом малым. Аналогия, думаю, ясна. :)

Предположим программа с dt_min попала в условия в которых программа с остановом остановилась бы. Что будет?
Программа с dt_min пойдёт дальше. И выведёт пользователю числа не имеющие никакого отношения к реальности. Разница будет не метра или два. Высота может стать отрицательной, например, или 10^100. А главное, что программа сделает это молча. Пользователь потеряет веру в её точность. А эта вера дорого стоит, её надо беречь.

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

Лунолёты 80-х считали с ∆t = 22 секунды по Эйлеру, не краснея. Числа были не совсем точными, но имели отношения к реальности. Наша программа, идя с ∆t = dt_min и обладая методом 4-го порядка, не будет завышать результаты в разы. Особенно при работе с системой, где нет ничего, сложнее возведения в квадрат и синуса гладко меняющегося (а то и постоянного) угла.

Когда будем шлифовать программу, добавлю проверку с аварийным остановом вместо Kmax, чтобы пользователь сам решал, идти ему выть на Луну от бессилия — или нажать С/П, получив (как раньше) осмысленный результат, но без гарантий точности. Сейчас этот блок лишний, он усложнит и без того не самую очевидную программу. У нас этих проверок будет ещё много — посадка, перегрузка. Может быть, даже стыковка.

Пока алгоритм не отлажен, веру цифрам всё равно культивировать нельзя. Каждый результат проверяем и перепроверяем.

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

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

После исправления бага Рунге-Кутта летает! И как летает!

точность e 1e-6
Полный бак m 3500
Манёвр 84 2100 74

Потребовался 41 шаг! Не один шаг кроме самого первого не был отменён. Размер шага в процессе манёвра уменьшался с двух секунд до полутора.
Экстраполяция делалась 123 раза.

точность e 1e-9
Шаг менялся от трети до четверти секунды.

Для справки: такая ошибка гарантирует, что через 10 суток полёта вы будете знать свои координаты с точностью до миллиметра.

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

Ещё тест:

Стартовал с полным баком, вышел на орбиту, сделал виток:

m 3500

84 2100 74
0 0 2086
90 47 2
0 0 7318

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

Ведь это означает, что Дорманд-Принс можно будет использовать для полётов реального времени — и точностью не придётся жертвовать!
Моя одиночная разработка, где я брал ∆t из таймера «Электроники МК-152», использовала плоскую декартовую систему координат и метод Эйлера.

> Например, стартуете с Байконура, как планируете попасть в плоскость орбиты Луны?

Я не пилот, и с реальными полётами знаком мало. Но Байконур раз в сутки максимально приближается к плоскости эклиптики (был бы южнее тропика, то он её бы пересекал два раза в сутки). Стартуем в этот момент по направлению вращения Земли. Полученная траектория будет иметь минимальный угол отклонения от плоскости эклиптики. Далее её можно будет скорректировать в плоскость эклиптики (в узлах где две плоскости пересекаются), а потом отправляться на Луну.

Это всё мои фантазии. Рассказал бы кто про реальность.

При обратном полёте всё проще, если нам безразлично, куда садится. Приземляемся куда нибудь между тропиками и ладно.

Извиняюсь, пропустил ваше сообщение, только сейчас заметил.
Тут важен наклон плоскости орбиты Луны относительно экватора. Он всегда меньше 28.6 градусов (23.5 - наклон экватора к эклиптике, 5.1 градусов - наклон орбиты Луны к эклиптике). Широта Байконура - 45 градусов, поэтому минимальный наклон орбиты, проходящей через эту точку, к экватору - 45 градусов. Стартовав с Байконура, чтоб выйти в плоскость Луны придется потратить дополнительно уйму топлива на изменение плоскости орбиты. Никто так не делает. Советские КА всегда летали вне плоскости орбиты Луны.

Модель с прямоугольными координатами позволяет добавить третье измерение (будет 6 координат в векторе). Все измения коснутся только функции g и расчёта ошибки. Надо придумать ввод и летайте на здоровье. Туда же можно добавить несколько притягивающих тел (только функция g). Можно будет повторить лунные миссии России и США.

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

Ну, система координат может быть и сферической...

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

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

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

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

Комментарий был перемещен и теперь находится здесь.

Отличная книга! По мере сил решения выписываю тут:
http://young-astronaut.livejournal.com/1062.html

---------------------------
Истина где-то рядом
www.litres.ru/vitaliy-samurov/dozvonitsya-do-devy/

Набрел на пресс-конференцию, которую проводила фирма Princeton Satellite Systems.
Они замахнулись на серьезный шаг: построить космический аппарат, оснащенный прямоточным термоядерным реактором-двигателем их разработки. Цель - доставить 1000 кг полезной нагрузки на орбиту Плутона меньше, чем за 4 года, и посадить спускаемый аппарат на поверхность Плутона.
(для сравнения - "Новые Горизонты" осуществил пролет возле Плутона через 9 лет после запуска. Также надо помнить, что пролет и выход на орбиту - совершенно разные вещи по требуемой энергетике).

Энерговооруженность двигателя - 2 МВатта электроэнергии во время всей миссии. Поэтому спускаемый аппарат будет получать около 30 КВатт энергии с орбитера посредством лазера. Все это позволит транслировать на Землю видео в высоком разрешении, в реальном времени.

Заявка очень сильная. Хотелось бы верить, но очень сложно поверить.

Еще на их сайте я нашел приложения для айфона.
Среди них - игра Outer Space Defense:

Outer Space Defense is a game designed to teach relative orbit dynamics. This game poses a test of skill pitting your mastery of spacecraft orbital control and prowess with a laser against multiple attacking enemy satellites. See how long you last against a relentless enemy that launches satellites designed to take you out at a faster and faster pace!

Подробное описание математики в игре тут.

P.S. С интересом обнаружил, что иногда подседаю на другой орбитальный "Space Defender": игра "Буран" на эмуляторе ПК "Специалист" :)

---------------------------
Истина где-то рядом
www.litres.ru/vitaliy-samurov/dozvonitsya-do-devy/

Страницы