Численные методы

Численный FORTH

03.04.2021 16:09:23 |

Автор: admin

Первое впечатление

Форт и сейчас известен, главным образом, среди разработки встроенных систем, как что-то вроде необычайного высокоуровневого ассемблера, например, для микроконтроллеров - AmForth и Mecrisp. Однако, когда-то давным давно был известен в другой ипостаси - как язык программирования научных приложений.

Книги о Форте

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

Таунсенд К., Фохт Д. ПРОЕКТИРОВАНИЕ И ПРОГРАММНАЯ РЕАЛИЗАЦИЯ ЭКСПЕРТНХ СИСТЕМ НА ПЕРСОНАЛЬНХ ЭВМ. М.: Финансы и статистика, 1990.

Я прочёл это и был впечатлён. Вот три хорошо знакомые мне книги:

Языки программирования в книжках, соответственно - Бейсик, Фортран и Форт! В книге Т. Тоффоли:

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

Физически CAM-6 состоит из модуля, который вставляется в один разъем IBM-PC (XT, AT или совместимых с ними моделей), и управляющего программного обеспечения, работающего в среде PC-DOS2. В то время как этот легко доступный головной компьютер обеспечивает размещение, экранирование, электропитание, дисковую память, монитор и стандартную операционную среду, вся действительная работа по моделированию клеточных автоматов с очень высокой скоростью совершается самим модулем с быстродействием, сравнимым (для этого частного приложения) с быстродействием CRAY-1. Управляющее программное обеспечение для CAM-6 написано на FORTH и работает на IBM-PC с памятью 256 К. Это программное обеспечение дополнено рядом готовых приложений и демонстрационных примеров и включает полный аннотированный список источников.

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

Очень интересно! Убедили, я тоже хочу попробовать Форт! Ощутить, как это, писать научное приложение на уровне 80-х или ранних 90-х. Говорят, что каждый фортер пишет свой Форт, но вряд ли это кому-то еще интересно, так что, пожалуй, попробую воспользоваться каким-то готовым Фортом и написать, скажем, программу для численного испытания... А, пускай будет это: C. Clay Marston and Gabriel G. BalintKurti. The Fourier grid Hamiltonian method for bound state eigenvalues and eigenfunctions // J. Chem. Phys. 91, 3571 (1989); doi: 10.1063/1.456888

Как будто на дворе 1989 год, компилятор Фортрана где-то на большой машине в институте, а у меня просто нечто под руками, никаких там Matlab в помине нет. Компьютер - чуть больше по своим возможностям, чем калькулятор.

Знакомство

Кстати, а вот и он - HP 35s. Программы в нём - шитый код (ниже рисунок из Руководства к калькулятору). Перейти с калькулятора на Форт очень просто - Форт тоже использует тот самый шитый код.

Калькулятор HP 35s. Шитый код.

В самом деле, каждое слово (термин Форта, аналог процедур или функций) состоит из других слов и так до самого нижнего уровня, где слова-примитивы это просто машинный код.

see normcdf: NORMCDF       flit 1.41421 F/ ERF F1.0 F+ flit .500000 F* ;oksee erf: ERF           FDUP FSIGN FSWAP FDUP F* FDUP flit .147000 F*                FDUP flit 1.27324 F+ FSWAP F1.0 F+ F/ F*                FNEGATE FEXP FNEGATE F1.0 F+ FSQRT F* ;oksee fsign: FSIGN         F0< DUP INVERT - S>F ;oksee dupDUP IS CODE    ( $4012D8 53 )                push    ebx

В Форте всё слова. Константа - слово, что оставляет на стеке число. Переменная оставляет на стеке адрес, где лежит её значение. Есть много разных книг и описаний Форта, не буду их повторять, однако, считаю, что одно из лучших описаний того, что и как происходит под капотом дано здесь - http://rigidus.ru/

По моему впечатлению, Форт язык прагматичный, а не парадигматичный. Это постфиксный язык? С одной стороны, да, с другой - не совсем. Вот скажем, если последовательно придерживаться постфиксной записи, то определения должны были бы быть наподобие следующего:(word) cvn { moveto show } def или { moveto show } /S exch def где ключевое слово def (определить) стоит после самого определения. Так делается в Postscript, но не в Форт. В Форте было бы : word moveto show ; - за двоеточием следует определяемый термин, точка с запятой завершает определение. Почему так? А так проще. Двоеточие переводит интерпретатор Форт в состояние компиляции STATE=-1 (true в Форте), и текст определения считывается слева на право, точка с запятой в состояние интерпретации (выполнения) STATE=0.

Многие вещи в Форте сделаны с упором на простоту реализации, а не стандартизации. Чак Мур, создатель Форта, плохо относится к стандартизации. Его принцип - тебе нужно, ты и сделай. Не спекулируй, не оставляй зарубки для удобства будущего расширения функционала. Сейчас удобно и просто делать ТАК, вот ТАК и делай САМ И ПРЯМО СЕЙЧАС. Короче говоря, Форт это не совсем язык, это метод решения конкретных проблем. Вот его кредо:

Keep it simple
Do not speculate
Do it yourself ~~motherf*cker~~

Не ждите от Форта серебряной пули, крутого фреймворка, полной библиотеки и продуманного кем-то API. Вот вам мнения, целый зоопарк реализаций и рецептов - да поможет Вам здравый смысл и базовые три принципа.

Стеквилибристика

Положим, что мы захотели реализовать вычисление гамма-функции (точнее, её логарифма) на Форте. Более того, у нас имеется сопроцессор типа Intel 8087 - у него стековая архитектура, очень кстати для Форта! Воспользуемся приближением Ланцоша и запишем:

: LNGAMMA ( x -- ln(Г(x) )\ Takes x > 0.0 and returns natural logarithm of Г(x).  FDUP  3.0E F+   2.23931796330267E FSWAP F/ FOVER 2.0E F+ -27.0638924937115E  FSWAP F/ F+ FOVER 1.0E F+  41.4174045302371E  FSWAP F/ F+     2.5066284643656E F+ FLN FSWAP FDUP  4.15E F+ FLN FOVER 0.50E F+ F* FOVER 4.15E F+ F-         FROT F+ FSWAP FLN F- ;

Да, выглядит не очень читаемо - не все формулы хорошо укладываются в операции со стеком и чем сложнее - тем больше будет слов FDUP, FROT, FOVER... пока не настанет ситуация, когда на стеке 4 и более чисел. Тогда всё, приплыли. Печальная история о том, как такое случается изложена в одном блоге.

Обычный выход из этого положения - локальные переменные. Да, это Форт, их можно реализовать разными способами. Например, так : lngamma { f: x } это способ gforth. Или так : lngamma {: f: x :} способ VFX Forth. Локальные переменные достаточно сложная концепция, включающая в себя область видимости, время жизни и прочее. Вот уже придется нарушить первый принцип Форт и превратить его в Си-подобный язык?

Гиперстатическое глобальное окружение

Строго говоря, глобальные переменные Форта это совсем не тоже самое, что в других языках. Чтобы было понятнее, вот пример рабочего кода:

variable apples  ok: +apples apples +! ;  ok: apples ." You have " apples @ . ." apples." cr ; okapples You have 0 apples. ok5 +apples  okapples You have 5 apples. ok

Здесь переменная apples переопределена словом apples, которое сообщает текущее количество яблок. Однако, слово +apples работает как положено, увеличивая счётчик. Внутри слова +apples ссылка на адрес счётчика, а не имя переменной. Так мы можем изменить любое определение не затрагивая работу ранних определений. Например, нам нужна переменная X. Определим:

variable &x: x &x @ ;: (x) &x ! ;: cube (x)  x x x * * ;variable &x: x &x @ ;: (x) &x ! ;: square (x)  x x * ;  3 square . 9  ok3 cube . 27  ok

Слова cube и square работают как положено. Повторяющиеся слова &x, x, (x) можно спрятать за определяющим словом, наподобие того, как предложено тут, см. пост FORTH: Самоопределяющиеся слова.

Определение группы переменных в F-PC Forth 3.60

FLOAD FFLOAT.SEQFLOAD EVAL.SEQ: COMPARE ( c-addr1 u1 c-addr2 u2 -- n )   ROT   2DUP U< IF DROP COMPARE DUP 0= IF DROP  1 THEN  EXIT THEN   2DUP U> IF NIP  COMPARE DUP 0= IF DROP -1 THEN  EXIT THEN   DROP COMPARE ;: REFILL ( -- f )  \ CORE version for user input device and string only   loading @      IF  ( file )                    false EXIT THEN   'tib @ sp0 @ = IF  ( user input device ) query  true EXIT THEN   ( EVALUATE )  false ;MACRO: ++ PAD +PLACE ;: (VARIABLE)  " VARIABLE &" PAD PLACE 2DUP ++  "  : (" ++ 2DUP ++ " ) &" ++   2DUP   ++ "  ! ;" ++  "  : "  ++ 2DUP ++ "   &" ++ ( NAME ) ++ "  @ ;" ++  PAD COUNT EVAL ;: (FVARIABLE)  " FVARIABLE &" PAD PLACE 2DUP ++  "  : (" ++ 2DUP ++ " ) &" ++   2DUP   ++ "  F! ;" ++  "  : "  ++ 2DUP ++ "   &" ++ ( NAME ) ++ "  F@ ;" ++  PAD COUNT EVAL ;: REFILL-AT-EOL? ( S: -- FLAG )  SOURCE NIP >IN @ > DUP 0= IF DROP REFILL THEN ;: VARIABLES(  BEGIN BL WORD COUNT 2DUP " )" COMPARE  WHILE REFILL-AT-EOL?  WHILE (VARIABLE)  REPEAT  THEN 2DROP ;: FVARIABLES(  BEGIN BL WORD COUNT 2DUP " )" COMPARE  WHILE REFILL-AT-EOL?  WHILE (FVARIABLE)  REPEAT  THEN 2DROP ;

Так оно будет использоваться так:

 \ Объявление переменных VARIABLES( MAXIT )FVARIABLES( ACCURACY UNLIKELY-VALUE )\ Присвоение значений-1.11E30 (UNLIKELY-VALUE) 1.0E-9  (ACCURACY)      50 (MAXIT)      \ Положить значение переменной на стекMAXIT . 50 ok

Польза от выбранной нотации, по сравнению с обычными переменными Форта, заключается в том, что гораздо чаще с переменной значение считывается, а не записывается. Таким образом, вместо x @ y @ + z ! будет x y + (z), и многочисленные @ и f@ будут факторизованы.

F-PC Forth

Для аутентичного погружения в старые добрые времена IBM PC AT, MS DOS и всё такое, я выбрал F-PC Forth. Скачать fpc36.zip, распаковать, запускать под dosbox. Работает всё из коробки, легко и просто.

F-PC Forth - стартовый экран.

Это полноценная IDE, с редактором кода, отладчиком и интерактивной справкой. Удобные IDE не только Borland делал.

больше ретро экранов F-PC Forth 3.60

Интерактивная справка.

Пошаговый отладчик кода. REPL

REPL

Написал в этом старом Форте код для поиска корней уравнения методом Риддера. Для сравнения здесь есть код на Julia.

Поиск корня уравнения методом Риддера на F-PC Forth 3.60

DEFER F(X) VARIABLES( MAXIT )FVARIABLES( XL XM XH XNEW  FL FM FH FNEW  S RESULT ACCURACY UNLIKELY-VALUE )-1.11E30 (UNLIKELY-VALUE) 1.0E-9  (ACCURACY)      50 (MAXIT): FSIGN ( R1 -- R1 ) F0< DUP NOT - IFLOAT ;: F~ ( R1 R2 R3 -- FLAG ) F-ROT F- FABS F> ;: ROOT-NOT-BRACKETED? ( FL FH -- FLAG )  FDUP   F0< FOVER  F0> AND  ( FB ) F0> ( FA ) F0< AND OR NOT ;: RIDDER ( R1 R2 -- R1 ) (XH) (XL)  XL F(X) (FL)  XH F(X) (FH)  FL F0= IF XL EXIT THEN  FH F0= IF XH EXIT THEN  FL FH ROOT-NOT-BRACKETED?  IF ABORT" ROOT MUST BE BRACKETED IN ZRIDDR" THEN  UNLIKELY-VALUE (RESULT) FALSE  MAXIT 0  DO    XL XH F+ 2.0E F/ (XM)  XM F(X) (FM)    FM FDUP F* FL FH F* F- FSQRT (S)    S F0=    IF RESULT TRUE LEAVE THEN    FL FH F- FSIGN XM XL F- F* FM F* S F/ XM F+ (XNEW)    XNEW RESULT ACCURACY F~    IF RESULT TRUE LEAVE THEN    XNEW (RESULT)  XNEW F(X) (FNEW)    FNEW F0=    IF RESULT TRUE LEAVE THEN    FNEW FSIGN FM F* FM F= NOT    IF XM (XL)  FM (FL)  RESULT (XH)  FNEW (FH)    ELSE FNEW FSIGN FL F* FL F= NOT      IF RESULT (XH)  FNEW (FH) THEN      FNEW FSIGN FH F* FH F= NOT      IF RESULT (XL)  FNEW (FL) THEN    THEN    XL XH ACCURACY F~    IF RESULT TRUE LEAVE THEN  LOOP  IF RESULT DROP  ELSE ." ZRIDDR EXCEED MAXIMUM ITERATIONS" DROP THEN ;: FUNC FDUP FEXP FSWAP -5.0E F* 3.0E F+ F+ ;' FUNC IS F(X)1.25E 1.6E RIDDER F.

Кажется, получилось вполне читаемо, особенно если сравнивать с языками тех времен: BASIC, Fortran 77, Pascal.

Структуры, массивы, матрицы

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

\ Structures: structure:  create ( structure name ) here 0 0 ,  does> @ ;: +field  create ( field name ) over , +  does> @ + ;: (cell)   aligned 1 cells  +field ;: (float) faligned 1 floats +field ;: end-structure ( addr size -- ) swap ! ;

Здесь я уже переключился на более современный Форт, стандарта 1994. В сущности, F-PC может быть дополнен до этого стандарта, и код ANS Forth 94 поддерживается современными компиляторами, например, win32forth, Gforth. Следуя духу ретро, я писал код в win32forth.

Win32forth IDE

У него тоже есть IDE и другие удобства, работает под Windows (под wine запускается без проблем). Структуры полезны при определении векторов и матриц, например:

\ Arraysstructure: array-structure  (cell) .array-data  (cell) .array-type  (cell) .array-sizeend-structure: array: ( size -- )  create  0 here .array-data ! here .array-type ! here .array-size !  array-structure allot ;: array-allocate ( vec -- )  >r r@ .array-size @ r@ .array-type @ * allocate throw r> .array-data ! ;: array-free ( vec -- )  >r r@ .array-data @ free throw 0 r> .array-data ! ;: array-element ( i vec -- *vec[i] )  >r r@ .array-type @ * r> .array-data @ + ;

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

code fs-array-element  pop eax  mov ebx, [ebx]  lea ebx, [ebx] [eax*8]next c;

Существует библиотека математических функций Форта - The Forth Scientific Library Project, впрочем, там всё равно нет реализации алгоритма вычисления собственных значений симметричной действительной матрицы. Do it yourself! Берем книгу Голуб, Ч. Ван Лоун. Матричные вычисления и реализуем алгоритм Якоби.

Циклический метод Якоби, псевдокод.

\ Cyclic Jacobi. Algorithm 8.5.3\ Golub & Van Loan, Matrix Computationsfvariables( cos sin EPS ) variables( M EV MAXROT )1.0e-10 (EPS)     50 (MAXROT): eig! (EV) (M)  EV matrix-set-identity!  MAXROT 0  do    M off-diagonal-norm EPS f<    if unloop exit then    M .matrix-rows @ 0    do M .matrix-cols @ i 1+       ?do i j M sym.schur2 (sin) (cos)        cos sin i j  M jacobi.rot'        cos sin i j  M jacobi.rot        cos sin i j EV jacobi.rot      loop    loop  loop  ." jacobi not converged" ;

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

Финал

Настало время решать поставленную задачу, а именно, согласно статье C. Clay Marston and Gabriel G. BalintKurti. The Fourier grid Hamiltonian method for bound state eigenvalues and eigenfunctions // J. Chem. Phys. 91, 3571 (1989); doi: 10.1063/1.456888 реализовать метод и посчитать, допустим, уровни энергии осциллятора Морзе. Иными словами, превратить рассуждения из фрагмента статьи ниже в практику.

The Fourier grid Hamiltonian method for bound state eigenvalues and eigenfunctions.

Последовательность действий:

Вычисляем суммы из уравнения (26)
Формируем тёплицеву матрицу H
Создаем дискретную сетку по оси X
Табулируем потенциал V(x) в точках сетки
Добавляем табулированный потенциал к диагонали матрицы H
Находим собственные числа и собственные векторы (это и есть решение)

Суть метода Fourier grid Hamiltonian (FGH) умещается в два определения Форта:

\ Equation 26fvariables( l d/N ): sum (d/N)  1.0e (l)  0.0e ( N ) 1 rshift 0  do [ 2.0e fpi f* ] fliteral     l d/N  f*  f* fcos l f**2 f* f+     l 1.0e f+ (l)  loop ; variables( diags n )fvariables( dx 1/n ): FGH! (diags) (dx)  diags .array-size @ (n)  n s>f 1/f (1/n)  [ -8.0e fpi f**2 f* ] fliteral  1/n fdup fdup f* f* f* dx f**2 1/f f*  n 0 do i s>f 1/n f* n sum fover f* i diags fa!  loop fdrop ;

Дальше идёт обычный boilerplate code как объявления матриц, векторов, инициализация элементов, затем поиск собственных значений и распечатка результатов. На этом этапе стоит как следует поиграть с параметрами и воспользоваться интерактивностью Форта. Мы же как будто в прошлом, так? Никакого python/numpy, Matlab и Julia - просто усовершенствованный калькулятор.

Результат вычислений.

Кому интересно, смотрите, я выложил код на гитхабе.

Заключение

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

Это всё философия, конечно. Однако, я могу себе представить какой-то численный Форт и сейчас, в наше время. Он может оказаться где то глубоко в оборудовании хитрого спектрометра, детектора... Было бы интересно узнать, где.

Подробнее..

Категории: Ненормальное программирование , Квантовая механика , Численные методы , Forth

Перенос молекулярной динамики на CUDA. Часть II Суммирование по Эвальду

07.07.2020 20:09:59 |

Автор: admin

В предыдущей статье мы обсудили основу метода молекулярной динамики, в том числе вычисление энергии и сил взаимодействия между частицами с заданными парными потенциалами. А что, если частицы обладают некоторым электрическим зарядом? Например, в том случае, если мы моделируем кристалл поваренной соли, состоящий из ионов Na⁺ и Cl^-. Или водный раствор, содержащий те или иные ионы. В этом случае, кроме парных потенциалов типа Леннарда-Джонса между ионами действуют силы электростатического взаимодействия, т.е. закон Кулона. Энергия такого взаимодействия для пары частиц i-j равна:

$E=C\frac{q_iq_j}{r_{ij}},$

где q заряд частицы, r_ij расстояние между частицами, С некоторая постоянная, зависящая от выбора единиц измерения. В системе СИ это $\frac{1}{4\pi\epsilon_0}$ , в СГС 1, в моей программе (где энергия выражена в электронвольтах, расстояние в ангстремах, а заряд в элементарных зарядах) C примерно равно 14.3996.

Ну и что, скажете вы? Просто добавим соответствующее слагаемое в парный потенциал и готово. Однако, чаще всего в МД моделировании используют периодические граничные условия, т.е. моделируемая система со всех сторон окружена бесконечным количеством её виртуальных копий. В этом случае каждый виртуальный образ нашей системы будет взаимодействовать со всеми заряженными частицами внутри системы по закону Кулона. А поскольку Кулоновское взаимодействие убывает с расстоянием очень слабо (как 1/r), то отмахнуться от него так просто нельзя, сказав, что с такого-то расстояния мы его не вычисляем. Ряд вида 1/x расходится, т.е. его сумма, в принципе, может расти до бесконечности. И что же теперь, миску супа не солить? Убьёт электричеством?

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

$E_{real} = С\sum\limits^{N-1}_i{\sum\limits^N_{j=i+1}{ \frac{q_iq_j}{r_{ij}}\operatorname{erfc}(\alpha r_{ij}) }},$

$E_{rec} = С\frac{2\pi}{V}\sum_{\vec{k}\ne0}(\exp(-k^2/4\alpha^2)/k^2 \sum_jq_i\exp(i\vec{k}\vec{r}_j)),$

$E_{const} = -\frac{C}{V\alpha} \sum_iq_i^2,$

где параметр, регулирующий соотношение вычислений в прямом и обратном пространствах, k вектора в обратном пространстве по которым идёт суммирование, V объём системы (в прямом пространстве). Первая часть (E_real) является короткодействующей и вычисляется в том же цикле, что и другие парные потенциалы, смотри функцию real_ewald в предыдущей статье. Последний вклад (E_сonst) является поправкой на самовзаимодействие и часто называется постоянной частью, поскольку не зависит от координат частиц. Её вычисление тривиально, поэтому мы остановимся только на второй части Эвальдовой суммы (E_rec), суммировании в обратном пространстве. Естественно, во времена вывода Эвальда молекулярной динамики не было, кто впервые использовал этот метод в МД найти мне не удалось. Сейчас любая книга по МД содержит его изложение как некий золотой стандарт. К книге Аллена даже прилагается пример кода на фортране. К счастью, у меня остался код, написанный когда-то на С для последовательной версии, осталось только его распараллелить (я позволил себе опустить некоторые объявления переменных и другие несущественные детали):

void ewald_rec(){    int mmin = 0;    int nmin = 1;    // массивы где хранятся iexp(x[i] * kx[l]),    double** elc;    double** els;    //... iexp(y[i] * ky[m]) и    double** emc;    double** ems;    //... iexp(z[i] * kz[n]),    double** enc;    double** ens;    // временные массивы для произведений iexp(x*kx)*iexp(y*ky)    double* lmc;    double* lms;    // и для q[i] * iexp(x*kx)*iexp(y*ky)*iexp(z*kz)    double* ckc;    double* cks;    // ПРЕДВАРИТЕЛЬНОЕ ЗАПОЛНЕНИЕ МАССИВОВ    eng = 0.0;    for (i = 0; i < Nat; i++)   // цикл по атомам    {        // emc/s[i][0] и enc/s[i][0] уже заполнены на этапе инициализации        // в массив elc/s нужно обновить, см. далее        elc[i][0] = 1.0;        els[i][0] = 0.0;        // iexp(kr)        sincos(twopi * xs[i] * ra, els[i][1], elc[i][1]);        sincos(twopi * ys[i] * rb, ems[i][1], emc[i][1]);        sincos(twopi * zs[i] * rc, ens[i][1], enc[i][1]);    }    // заполняем следующие элементы массива emc/s[i][l] = iexp(y[i]*ky[l]) итеративно, используя комплексное умножение    for (l = 2; l < ky; l++)        for (i = 0; i < Nat; i++)        {            emc[i][l] = emc[i][l - 1] * emc[i][1] - ems[i][l - 1] * ems[i][1];            ems[i][l] = ems[i][l - 1] * emc[i][1] + emc[i][l - 1] * ems[i][1];        }    // заполняем следующие элементы массива enc/s[i][l] = iexp(z[i]*kz[l]) итеративно, используя комплексное умножение    for (l = 2; l < kz; l++)        for (i = 0; i < Nat; i++)        {            enc[i][l] = enc[i][l - 1] * enc[i][1] - ens[i][l - 1] * ens[i][1];            ens[i][l] = ens[i][l - 1] * enc[i][1] + enc[i][l - 1] * ens[i][1];        }    // ГЛАВНЙ ЦИКЛ ПО ВСЕМ K-ВЕКТОРАМ:    for (l = 0; l < kx; l++)    {        rkx = l * twopi * ra;         // записываем exp(ikx[l]) в ikx[0] для сохранения памяти        if (l == 1)            for (i = 0; i < Nat; i++)            {                elc[i][0] = elc[i][1];                els[i][0] = els[i][1];            }        else if (l > 1)            for (i = 0; i < Nat; i++)            {                // iexp(kx[0]) = iexp(kx[0]) * iexp(kx[1])                x = elc[i][0];                elc[i][0] = x * elc[i][1] - els[i][0] * els[i][1];                els[i][0] = els[i][0] * elc[i][1] + x * els[i][1];            }        for (m = mmin; m < ky; m++)        {            rky = m * twopi * rb;            // заполняем временный массив произведением iexp(kx*x[i]) * iexp(ky*y[i])            if (m >= 0)                for (i = 0; i < Nat; i++)                {                    lmc[i] = elc[i][0] * emc[i][m] - els[i][0] * ems[i][m];                    lms[i] = els[i][0] * emc[i][m] + ems[i][m] * elc[i][0];                }            else // для отрицательных значений m используем комплексное сопряжение:                for (i = 0; i < Nat; i++)                {                    lmc[i] = elc[i][0] * emc[i][-m] + els[i][0] * ems[i][-m];                    lms[i] = els[i][0] * emc[i][-m] - ems[i][-m] * elc[i][0];                }            for (n = nmin; n < kz; n++)            {                rkz = n * twopi * rc;                rk2 = rkx * rkx + rky * rky + rkz * rkz;                if (rk2 < rkcut2) // используем радиус обрезания                {                    // вычисляем сумму (q[i]*iexp(kr[k]*r[i])) - зарядовую плотность                    sumC = 0; sumS = 0;                    if (n >= 0)                        for (i = 0; i < Nat; i++)                        {                            //считываем заряд частицы                            ch = charges[types[i]].charge;                            ckc[i] = ch * (lmc[i] * enc[i][n] - lms[i] * ens[i][n]);                            cks[i] = ch * (lms[i] * enc[i][n] + lmc[i] * ens[i][n]);                            sumC += ckc[i];                            sumS += cks[i];                        }                    else // для отрицательных индексов используем комплексное сопряжение:                        for (i = 0; i < Nat; i++)                        {                            //считываем заряд частицы                            ch = charges[types[i]].charge;                            ckc[i] = ch * (lmc[i] * enc[i][-n] + lms[i] * ens[i][-n]);                            cks[i] = ch * (lms[i] * enc[i][-n] - lmc[i] * ens[i][-n]);                            sumC += ckc[i];                            sumS += cks[i];                        }                    //наконец вычисляем энергию и силы                    akk = exp(rk2 * elec->mr4a2) / rk2;                    eng += akk * (sumC * sumC + sumS * sumS);                    for (i = 0; i < Nat; i++)                    {                        x = akk * (cks[i] * sumC - ckc[i] * sumS) * C * twopi * 2 * rvol;                        fxs[i] += rkx * x;                        fys[i] += rky * x;                        fzs[i] += rkz * x;                    }                }            } // end n-loop (over kz-vectors)            nmin = 1 - kz;        } // end m-loop (over ky-vectors)        mmin = 1 - ky;    }  // end l-loop (over kx-vectors)   engElec2 = eng * С * twopi * rvol;}

Пара пояснений к коду: функция считает комплексную экспоненту (в комментариях к коду она обозначена iexp, чтобы убрать мнимую единицу из скобок) от векторного произведения k-вектора на радиус-вектор частицы для всех k-векторов и для всех частиц. Эта экспонента домножается на заряд частицы. Далее вычисляется сумма таких произведений по всем частицам (внутренняя сумма в формуле для E_rec), у Френкеля она называется зарядовой плотностью, а у Блинова структурным фактором. Ну а далее, на основании этих структурных факторов вычисляется энергия и силы, действующие на частицы. Компоненты k-векторов (2*l/a, 2*m/b, 2*n/c) характеризуются тройкой целых чисел l, m и n, которые и пробегаются в циклах до заданных пользователем пределов. Параметры a, b и c это размеры моделируемой системы в измерениях x, y и z соответственно (вывод верен для системы с геометрией прямоугольного параллелепипеда). В коде 1/a, 1/b и 1/с соответствуют переменным ra, rb и rc. Массивы под каждую величину представлены в двух экземплярах: под действительную и мнимую части. Каждый следующий k-вектор в одном измерении получается итеративно из предыдущего путем комплексного умножения предыдущего на единичный, чтобы каждый раз не считать синус с косинусом. Массивы emc/s и enc/s заполняются для всех m и n, соответственно, а массив elc/s значение для каждого l>1 помещает в нулевой индекс по l в целях экономии памяти.
В целях распараллеливания выгодно вывернуть порядок циклов так, чтобы внешний цикл пробегался по частицам. И тут мы видим проблему распараллелить эту функцию можно только до вычисления суммы по всем частицам (зарядовой плотности). Дальнейшие вычисления опираются на эту сумму, а она будет рассчитана только когда все потоки закончат работу, поэтому придётся разбить эту функцию на две. Первая вычисляет считает зарядовую плотность, а вторая энергию и силы. Замечу, что во второй функции снова потребуется величина q_iiexp(kr) для каждой частицы и для каждого k-вектора, вычисленная на предыдущем этапе. И тут есть два подхода: либо пересчитать её заново, либо запомнить. Первый вариант требует больше времени, второй больше памяти (количество частиц * количество k-векторов * sizeof(float2)). Я остановился на втором варианте:

__global__ void recip_ewald(int atPerBlock, int atPerThread, cudaMD* md)// calculate reciprocal part of Ewald summ// the first part : summ (qiexp(kr)) evaluation{    int i;      // for atom loop    int ik;     // index of k-vector    int l, m, n;    int mmin = 0;    int nmin = 1;    float tmp, ch;    float rkx, rky, rkz, rk2;   // component of rk-vectors    int nkx = md->nk.x;    int nky = md->nk.y;    int nkz = md->nk.z;        // arrays for keeping iexp(k*r) Re and Im part    float2 el[2];    float2 em[NKVEC_MX];    float2 en[NKVEC_MX];    float2 sums[NTOTKVEC];          // summ (q iexp (k*r)) for each k    extern __shared__ float2 sh_sums[];     // the same in shared memory    float2 lm;     // temp var for keeping el*em    float2 ck;     // temp var for keeping q * el * em * en (q iexp (kr))    // invert length of box cell    float ra = md->revLeng.x;    float rb = md->revLeng.y;    float rc = md->revLeng.z;    if (threadIdx.x == 0)        for (i = 0; i < md->nKvec; i++)            sh_sums[i] = make_float2(0.0f, 0.0f);    __syncthreads();    for (i = 0; i < md->nKvec; i++)        sums[i] = make_float2(0.0f, 0.0f);    int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;    int N = min(id0 + atPerThread, md->nAt);    ik = 0;    for (i = id0; i < N; i++)    {        // save charge        ch = md->specs[md->types[i]].charge;        el[0] = make_float2(1.0f, 0.0f);    // .x - real part (or cos) .y - imagine part (or sin)                em[0] = make_float2(1.0f, 0.0f);        en[0] = make_float2(1.0f, 0.0f);        // iexp (ikr)        sincos(d_2pi * md->xyz[i].x * ra, &(el[1].y), &(el[1].x));        sincos(d_2pi * md->xyz[i].y * rb, &(em[1].y), &(em[1].x));        sincos(d_2pi * md->xyz[i].z * rc, &(en[1].y), &(en[1].x));        // fil exp(iky) array by complex multiplication        for (l = 2; l < nky; l++)        {             em[l].x = em[l - 1].x * em[1].x - em[l - 1].y * em[1].y;             em[l].y = em[l - 1].y * em[1].x + em[l - 1].x * em[1].y;        }        // fil exp(ikz) array by complex multiplication        for (l = 2; l < nkz; l++)        {             en[l].x = en[l - 1].x * en[1].x - en[l - 1].y * en[1].y;             en[l].y = en[l - 1].y * en[1].x + en[l - 1].x * en[1].y;        }        // MAIN LOOP OVER K-VECTORS:        for (l = 0; l < nkx; l++)        {            rkx = l * d_2pi * ra;                         // move exp(ikx[l]) to ikx[0] for memory saving (ikx[i>1] are not used)            if (l == 1)                el[0] = el[1];            else if (l > 1)                {                    // exp(ikx[0]) = exp(ikx[0]) * exp(ikx[1])                    tmp = el[0].x;                    el[0].x = tmp * el[1].x - el[0].y * el[1].y;                    el[0].y = el[0].y * el[1].x + tmp * el[1].y;                }            //ky - loop:            for (m = mmin; m < nky; m++)            {                rky = m * d_2pi * rb;                //set temporary variable lm = e^ikx * e^iky                if (m >= 0)                {                        lm.x = el[0].x * em[m].x - el[0].y * em[m].y;                         lm.y = el[0].y * em[m].x + em[m].y * el[0].x;                }                else // for negative ky give complex adjustment to positive ky:                {                        lm.x = el[0].x * em[-m].x + el[0].y * em[-m].y;                        lm.y = el[0].y * em[-m].x - em[-m].x * el[0].x;                }                //kz - loop:                for (n = nmin; n < nkz; n++)                {                    rkz = n * d_2pi * rc;                    rk2 = rkx * rkx + rky * rky + rkz * rkz;                    if (rk2 < md->rKcut2) // cutoff                    {                        // calculate summ[q iexp(kr)]   (local part)                        if (n >= 0)                         {                                ck.x = ch * (lm.x * en[n].x - lm.y * en[n].y);                                ck.y = ch * (lm.y * en[n].x + lm.x * en[n].y);                         }                        else // for negative kz give complex adjustment to positive kz:                         {                                ck.x = ch * (lm.x * en[-n].x + lm.y * en[-n].y);                                ck.y = ch * (lm.y * en[-n].x - lm.x * en[-n].y);                        }                        sums[ik].x += ck.x;                        sums[ik].y += ck.y;                                                // save qiexp(kr) for each k for each atom:                        md->qiexp[i][ik] = ck;                        ik++;                    }                } // end n-loop (over kz-vectors)                nmin = 1 - nkz;            } // end m-loop (over ky-vectors)            mmin = 1 - nky;        }  // end l-loop (over kx-vectors)    } // end loop by atoms    // save sum into shared memory    for (i = 0; i < md->nKvec; i++)    {        atomicAdd(&(sh_sums[i].x), sums[i].x);        atomicAdd(&(sh_sums[i].y), sums[i].y);    }    __syncthreads();    //...and to global    int step = ceil((double)md->nKvec / (double)blockDim.x);    id0 = threadIdx.x * step;    N = min(id0 + step, md->nKvec);    for (i = id0; i < N; i++)    {        atomicAdd(&(md->qDens[i].x), sh_sums[i].x);        atomicAdd(&(md->qDens[i].y), sh_sums[i].y);    }}// end 'ewald_rec' function__global__ void ewald_force(int atPerBlock, int atPerThread, cudaMD* md)// calculate reciprocal part of Ewald summ// the second part : enegy and forces{    int i;      // for atom loop    int ik;     // index of k-vector    float tmp;    // accumulator for force components    float3 force;    // constant factors for energy and force    float eScale = md->ewEscale;    float fScale = md->ewFscale;    int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;    int N = min(id0 + atPerThread, md->nAt);    for (i = id0; i < N; i++)    {        force = make_float3(0.0f, 0.0f, 0.0f);        // summ by k-vectors        for (ik = 0; ik < md->nKvec; ik++)        {            tmp = fScale * md->exprk2[ik] * (md->qiexp[i][ik].y * md->qDens[ik].x - md->qiexp[i][ik].x * md->qDens[ik].y);            force.x += tmp * md->rk[ik].x;            force.y += tmp * md->rk[ik].y;            force.z += tmp * md->rk[ik].z;        }        md->frs[i].x += force.x;        md->frs[i].y += force.y;        md->frs[i].z += force.z;    } // end loop by atoms    // one block calculate energy    if (blockIdx.x == 0)        if (threadIdx.x == 0)        {            for (ik = 0; ik < md->nKvec; ik++)                md->engCoul2 += eScale * md->exprk2[ik] * (md->qDens[ik].x * md->qDens[ik].x + md->qDens[ik].y * md->qDens[ik].y);        }}// end 'ewald_force' function

Надеюсь, вы мне простите, что я оставил комментарии на английском, код практически повторяет последовательную версию. Код даже стал читабельнее за счет того, что массивы потеряли одно измерение: elc/s[i][l], emc/s[i][m] и enc/s[i][n] превратились в одномерные el, em и en, массивы lmc/s и ckc/s в переменные lm и ck (пропала мерность по частицам, поскольку отпала необходимость хранить это для каждой частицы, промежуточный результат накапливается в shared memory). К сожалению, тут же возникла и проблема: массивы em и en пришлось задать статическими, чтобы не использовать глобальную память и не выделять память динамически каждый раз. Количество элементов в них определяется директивой NKVEC_MX (максимальное количество k-векторов по одному измерению) на этапе компиляции, а runtime используются только первые nky/z элементов. Кроме того, появился сквозной индекс по всем k-векторам и аналогичная директива, ограничивающая общее количество этих векторов NTOTKVEC. Проблема возникнет, если пользователю понадобится больше k-векторов, чем определено директивами. Для вычисления энергии предусмотрен блок с нулевым индексом, поскольку неважно какой именно блок выполнит этот расчет и в каком порядке. Тут может быть надо было использовать готовые функции FFT, раз уж метод основан на нём, но я так и не сообразил, как это сделать.
Ну а теперь попробуем что-нибудь посчитать, да тот же содиум хлорайд. Возьмём 2 тысячи ионов натрия и столько же хлора. Заряды зададим целыми, а парные потенциалы возьмём, например, из этой работы. Стартовую конфигурацию зададим случайно и слегка перемешаем её, рисунок 2а. Объём системы выберем так, чтобы он соответствовал плотности поваренной соли при комнатной температуре (2,165 г/см³). Запустим все это на небольшое время (10000 шагов по 5 фемтосекунд) с наивным учетом электростатики по закону Кулона и используя суммирование по Эвальду. Результирующие конфигурации приведены на рисунках 2б и 2в, соответственно. Вроде бы в случае с Эвальдом система чуть больше упорядочилась чем без него. Важно также, что флуктуации полной энергии с применением суммирования существенно уменьшились.

Рисунок 2. Начальная конфигурация системы NaCl (a) и после 10000 шагов интегрирования: наивным способом (б) и со схемой Эвальда (в).

Вместо заключения

Замечу, что структура, получаемая на рисунке, не соответствует кристаллической решетки NaCl, а скорее решетке ZnS, но это уже претензия к парным потенциалам. Учет же электростатики очень важен для молекулярно-динамического моделирования. Считается, что именно электростатическое взаимодействие ответственно за образование кристаллических решёток, поскольку действует на больших расстояниях. Правда с этой позиции сложно объяснить как при охлаждении кристаллизуются такие вещества как аргон.
Кроме упомянутого метода Эвальда, есть ещё и другие способы учета электростатики, смотрите, например, этот обзор.

Подробнее..

Категории: Научно-популярное , Параллельное программирование , Физика , Математическое моделирование , Химия , Cuda , Молекулярная динамика , Электростатика , Суммирование по эвальду , Численные методы

Равномерное перемещение объекта вдоль кривой

27.07.2020 18:22:43 |

Автор: admin

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

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

В нашем случае сплайн это функция, отображающая набор входных параметров (контрольные точки) и значение аргумента $inline$ (обычно принимающего значения от 0 до 1) в точку на плоскости или в пространстве. Получившаяся кривая это множество значений функции-сплайна при $t\in[0,1]$ .

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

Кубическая кривая Безье

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

Анимация отображения параметра t в точку кривой

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

Кусочный сплайн

С заданием и параметризацией маршрута разобрались, теперь переходим к основному вопросу. Чтобы наш условный самолёт двигался вдоль маршрута с постоянной скоростью, нам нужно в любой момент времени уметь вычислять точку на кривой в зависимости от пройденного вдоль этой кривой расстояния $inline$ , при этом имея лишь возможность вычислить положение точки на кривой по значению параметра $inline$ (функцию $inline$ ). Именно на этом этапе начинаются сложности.

Первой мыслью было сделать линейное отображение $len\in[0, Length] \Rightarrow t\in[0,1]$ и вычислить $inline$ от получившегося значения легко, вычислительно дёшево, в общем, то, что нужно.

Проблема данного способа очевидна сразу же на самом деле пройденная дистанция $inline$ зависит от $inline$ нелинейно, и чтобы в этом убедиться достаточно расставить по маршруту равномерно распределенные по значению $inline$ точки:

Равномерно распределенные по пути точки

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

Визуализация неправильной параметризации кривой

Обратившись за советом в поисковик, на stackoverflow и youtube, я обнаружил второй способ вычисления $inline$ , а именно представление кривой в виде кусочно-линейной функции (расчета набора равноудаленных друг от друга по кривой точек):

Представление кривой в виде кусочно-линейного сплайна

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

Пример кода

Источник

public Vector2[] CalculateEvenlySpacedPoints(float spacing, float resolution = 1)  {    List<Vector2> evenlySpacedPoints = new List<Vector2>();    evenlySpacedPoints.Add(points[0]);    Vector2 previousPoint = points[0];    float dstSinceLastEvenPoint = 0;    for (int segmentIndex = 0; segmentIndex < NumSegments; segmentIndex++)    {      Vector2[] p = GetPointsInSegment(segmentIndex);      float controlNetLength = Vector2.Distance(p[0], p[1]) + Vector2.Distance(p[1], p[2]) + Vector2.Distance(p[2], p[3]);      float estimatedCurveLength = Vector2.Distance(p[0], p[3]) + controlNetLength / 2f;      int divisions = Mathf.CeilToInt(estimatedCurveLength * resolution * 10);      float t = 0;      while (t <= 1)      {        t += 1f/divisions;        Vector2 pointOnCurve = Bezier.EvaluateCubic(p[0], p[1], p[2], p[3], t);        dstSinceLastEvenPoint += Vector2.Distance(previousPoint, pointOnCurve);        while (dstSinceLastEvenPoint >= spacing)        {          float overshootDst = dstSinceLastEvenPoint - spacing;          Vector2 newEvenlySpacedPoint = pointOnCurve + (previousPoint - pointOnCurve).normalized * overshootDst;          evenlySpacedPoints.Add(newEvenlySpacedPoint);          dstSinceLastEvenPoint = overshootDst;          previousPoint = newEvenlySpacedPoint;        }        previousPoint = pointOnCurve;      }    }    return evenlySpacedPoints.ToArray();  }

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

Вследствие чего я продолжил поиски и наткнулся на отличную статью Moving Along a Curve with Specified Speed, на основе которой строятся дальнейшие рассуждения.

Значение скорости можно вычислить как $\sigma(t) = |V(t)|=|\frac{dX}{dt}|$ , где $inline$ функция сплайна. Чтобы решить поставленную задачу, нужно найти функцию $inline$ , где $inline$ длина дуги от начальной точки до искомой, и это выражение устанавливает соотношение между $inline$ и $inline$ .

Используя замену переменной дифференцирования, мы можем записать

$\frac{dY}{dt}=\frac{dX}{ds}\frac{ds}{dt}.$

Учитывая, что скорость величина неотрицательная, получаем

$|\frac{dY}{dt}|=|\frac{dX}{ds}||\frac{ds}{dt}|=\frac{ds}{dt},$

так как $|\frac{dX}{ds}| = 1$ по условию неизменности модуля вектора скорости (вообще говоря, $|\frac{dX}{ds}|$ равен не единице, а константе, но для простоты выкладок примем эту константу равной единице).

Теперь получим соотношение между $inline$ и $inline$ в явном виде:

$s=g(t)=\int_{0}^{t}|\frac{dY(\tau)}{dt}|d\tau,$

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

Нас же интересует обратная операция, то есть вычисление значения аргумента $inline$ по заданной длине дуги $inline$ :

$t=g^{-1}(s).$

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

Метод образует последовательность приближений вида

$t_{i+1}=t_i-\frac{F(t_i)}{F'(t_i)},$

где

$F'(t)=\frac{dF}{dt}=\frac{dg}{dt}=|\frac{dY}{dt}|.$

Для вычисления $inline$ требуется вычислить $inline$ , вычисление которого, в свою очередь, требует численного интегрирования.

Выбор $t_0=\frac{s}{L}$ в качестве начального приближения теперь выглядит разумным (вспоминаем первый подход к решению задачи).

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

При использовании стандартного метода Ньютона потенциально может возникнуть проблема. Функция $inline$ неубывающая, так как её производная $inline$ неотрицательна. Если вторая производная $inline$ , то функцию называют выпуклой и метод Ньютона на ней гарантированно сходится к корню. Однако, в нашем случае, $inline$ может оказаться отрицательной, из-за чего метод Ньютона может сойтись за пределами диапазона определения $t\in[0,1]$ . Автор статьи предлагает использовать модификацию метода Ньютона, которая, в случае если очередной результат итерации методом Ньютона не попадает в текущий ограничивающий корень интервал, вместо него выбирает середину интервала (метод дихотомии). Вне зависимости от результата вычисления на текущей итерации, диапазон сужается, а значит, метод сходится к корню.

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

Код функции, вычисляющей t(s)

public float ArcLength(float t) => Integrate(x => TangentMagnitude(x), 0, t);private float Parameter(float length){  float t = 0 + length / ArcLength(1);  float lowerBound = 0f;   float upperBound = 1f;  for (int i = 0; i < 100; ++i)  {    float f = ArcLength(t) - length;    if (Mathf.Abs(f) < 0.01f)      break;    float derivative = TangentMagnitude(t);    float candidateT = t - f / derivative;    if (f > 0)    {      upperBound = t;      if (candidateT <= 0)        t = (upperBound + lowerBound) / 2;      else        t = candidateT;    }    else    {      lowerBound = t;      if (candidateT >= 1)        t = (upperBound + lowerBound) / 2;      else        t = candidateT;    }  }  return t;}

Код функции численного интегрирования

private static readonly (float, float)[] CubicQuadrature ={(-0.7745966F, 0.5555556F), (0, 0.8888889F), (0.7745966F, 0.5555556F)};public static float Integrate(Func<float, float> f, in float lowerBound, in float uppedBound){  var sum = 0f;  foreach (var (arg, weight) in CubicQuadrature)  {    var t = Mathf.Lerp(lowerBound, uppedBound, Mathf.InverseLerp(-1, 1, arg));    sum += weight * f(t);  }  return sum * (upperBound - lowerBound) / 2;}

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

Равномерно распределенные вдоль пути точки

Теперь самолет движется равномерно

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

I prefer the Newtons-method approach because generally t can be computed faster than with differential equation solvers.

В качестве заключения приведу таблицу времени расчета положения точки на представленной на скриншотах кривой в одном потоке на процессоре i5-9600K:

Кол-во вычислений	Среднее время, мс
100	0,62
1000	6,24
10000	69,01
100000	672,81

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

Подробнее..

Категории: Алгоритмы , C , Математика , Разработка игр , Геймдев , Unity , Численные методы , Кривые , Метод ньютона , Метод гаусса

	Русский
	English

Численные методы

Численный FORTH

Первое впечатление

Знакомство

Стеквилибристика

Гиперстатическое глобальное окружение

F-PC Forth

Структуры, массивы, матрицы

Финал

Заключение

Перенос молекулярной динамики на CUDA. Часть II Суммирование по Эвальду

Вместо заключения

Равномерное перемещение объекта вдоль кривой

Категории

Последние комментарии