Использование Джулии для решения высокоточных численных задач.

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

Доктор Трефетен не ожидал, что его вызов вызовет волну интереса. Более ста математиков со всего мира приняли участие в конкурсе и представили решения, причем 20 из них представили решения с высшим баллом в 100 правильных цифр. Чтобы узнать больше, введите в поисковой системе «SIAM 100-значный вызов», чтобы получить дополнительную информацию о конкурсе.

Хотя соревнование давно закончилось, решение любой из этих задач улучшит ваши навыки численного анализа. Некоторые из задач требуют продвинутых математических понятий, но большинство из них доступно любому, кто изучал математический анализ на уровне колледжа. Задача 2 хороша для начала, так как она относительно проста и не требует специальных навыков высокоточной арифметики. Использование типа BigFloat в Julia позволяет найти простое решение. Вот проблема:

A photon moving at speed 1 in the x-y plane starts at t = 0 at
    (x, y) = (1/2, 1/10)
heading due east. Around every integer lattice point (i, j) in the plane, a
circular mirror of radius 1/3 has been erected. How far from the origin is
the photon at t = 10?

Вам нужно немного узнать о том, как световой луч отражается от круглого зеркала, чтобы решить эту проблему. На приведенной ниже диаграмме луч света входит сверху слева и падает на зеркало в точке P, а затем выходит вверх и вправо. Пунктирная линия, обозначенная N, является линией нормали к поверхности в точке P. В данном случае это линия, проходящая через центр окружности и точку P. Если угол, который падающий луч образует с нормалью, равен θ, то угол исходящего луча также θ.

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

Самая интересная часть проблемы — определить, какое зеркало будет следующим, в которое ударит фотон. Метод грубой силы состоял бы в том, чтобы по очереди смотреть на каждое зеркало поблизости в поисках пересечений, но это неэффективно и не нужно. Вместо этого рассмотрим компоненты x и y скорости фотона, vx и vy, когда он покидает предыдущее зеркало. Мы соберем следующую информацию:

  • Is vx > 0?
  • Is vy > 0?
  • Is vx > vy?

Для конкретного примера выберем vx=0,5 и vy=0,866. Поскольку vx > 0, мы можем сразу исключить любое зеркало слева от текущего зеркала. Точно так же, поскольку vy › 0, мы можем исключить любое зеркало ниже текущего зеркала. На следующей диаграмме показана ситуация.

Красный круг — это место, где берет начало луч. Просто взглянув на знак компонентов скорости, мы убрали все серые зеркала.

Далее, поскольку vy › vx, фотон движется по относительно крутому пути, больше вверх, чем вправо. Это устраняет зеркала справа от текущего зеркала. Из 8 зеркал непосредственно вокруг текущих зеркал мы исключили 6 из них, оставив только зеркало сверху и зеркало сверху и справа. Луч может попасть в одну из них или пройти между ними. Если он проходит между ними, вы, вероятно, видите, что он должен в конечном итоге попасть в одно из зеркал в столбце справа от текущего зеркала. Следующая диаграмма показывает это.

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

Для решения этой проблемы я написал код Julia в блокноте Jupyter, который можно найти по адресу https://github.com/ciric50/julia-examples.git. Не стесняйтесь загружать его, но я настоятельно рекомендую вам сначала попробовать собственное решение. Так вы узнаете больше и, возможно, найдете алгоритм лучше, чем тот, который открыл я! Если вы хотите проверить свое решение на правильность, в моей записной книжке есть конечное число с гораздо большим количеством цифр, чем требуемые 10.