Квантовая физика на компьютере. Часть 2. Квантовая философия

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

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

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

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

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

Для этого мы используем метод расчёта матриц перехода, который называется flow theory. А матрицы перехода позволят построить “траекторию” частицы.

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

Напомню схему эксперимента:

Простой эксперимент с интерференцией. Фотоны вылетают из источника S, проходят через первое полупрозрачное зеркало, отражаются от промежуточных зеркал L, R, попадают на второе полупрозрачное зеркало, и, наконец, в один из двух детекторов T, U.

Исходными данными для расчёта служит начальное квантовое состояние $\psi^{(0)}$ и набор унитарных матриц $U^{(t)}$ , которые задают динамику квантового состояния:

$\psi^{(t)} = U^{(t)} \psi^{(t-1)}.$

Наше квантовое состояние очень простое — это две амплитуды:

$\psi^{(t)} = \begin{pmatrix} C_1 \\ C_2 \end{pmatrix},$

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

Наша цель — превратить $U^{(t)}$ и $\psi^{(t)}$ в матрицы перехода $\Gamma^{(t)}$ такие, что если на -м шаге вероятность пребывания в состоянии равна $P^{(t)}_1$ , а в состоянии 2 — $P^{(t)}_2$ , то

$\begin{pmatrix} P^{(t)}_1 \\ P^{(t)}_2 \end{pmatrix} = \Gamma^{(t)} \begin{pmatrix} P^{(t-1)}_1 \\ P^{(t-1)}_2 \end{pmatrix}.$

Как же их вычислить?

Согласно flow theory для этого надо рассчитать максимальный поток вот в таком направленном графе:

(На рёбрах написана их максимальная пропускная способность.) Тогда отношение потока в средних рёбрах (обозначенных $|U^{(t)}_{\cdot\cdot}|$ ) к соответствующим входным потокам ( $|\psi^{(t-1)}_{\cdot}|^2$ ) и даст элементы матрицы $\Gamma^{(t)}$ .

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

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

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

function build_flow_graph(start, finish, through) result(graph)  use, intrinsic :: iso_fortran_env, only: dp => real64  real(dp), intent(in), dimension(:) :: start  real(dp), intent(in), dimension(:) :: finish  real(dp), intent(in), dimension(:, :) :: through  real(dp), dimension(size(start) + size(finish) + 2, size(start) + size(finish) + 2) :: graph  integer :: i  integer :: j  if (size(start) /= size(through, 2)) error stop "hidden_flow:build_flow_graph:arg"  if (size(finish) /= size(through, 1)) error stop "hidden_flow_build_flow_graph:arg"  graph = 0.0_dp  do i = 1, size(start)    graph(2 + i, 1) = start(i)  end do  do j = 1, size(finish)    graph(2, 2 + size(start) + j) = finish(j)  end do  do i = 1, size(start)    do j = 1, size(finish)      graph(2 + size(start) + j, 2 + i) = through(j, i)    end do  end doend function

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

function find_shortest_path(graph, flow) result(path)  use, intrinsic :: iso_fortran_env, only: dp => real64  real(dp), intent(in), dimension(:, :) :: graph  real(dp), intent(in), dimension(:, :) :: flow  integer, dimension(:), allocatable :: path  <...>  if (size(graph, 1) /= size(graph, 2)) error stop "hidden_flow:find_shortest_path:arg"  if (size(flow, 1) /= size(flow, 2)) error stop "hidden_flow:find_shortest_path:arg"  if (size(flow, 1) /= size(graph, 1)) error stop "hidden_flow:find_shortest_path:arg"  if (size(graph, 1) < 2) error stop "hidden_flow:find_shortest_path:arg"

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

real(dp), parameter :: eps = 1.0e-10_dp

Главная идея алгоритма — использовать очередь (FIFO) узлов. Для очереди сделаем кольцевой буфер:

integer, dimension(size(graph, 1) + 1) :: queueinteger :: queue_head, queue_tail

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

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

logical, dimension(size(graph, 1)) :: visitedinteger, dimension(size(graph, 1)) :: prev

Для сборки пути неизвестной длины нужны переаллокации, значит нужно запоминать сколько памяти выделено и сколько из неё использовано, также понадобится второй временный буфер:

integer :: temp_capacityinteger :: temp_leninteger, dimension(:), allocatable :: temp

Общий алгоритм такой: 1) инициализируем очередь истоком, 2) пока очередь не исчерпается, достаём из неё узел, проходим по всем рёбрам из него в ещё не посещённые узлы и добавляем их в очередь, 3) когда дошли до стока, кратчайший путь найден, 4) если очередь исчерпалась, то пути нет.

queue(1) = 1queue_head = 1queue_tail = 2visited = .false.visited(1) = .true.prev(1) = 0q: do  if (queue_tail == queue_head) then    allocate(path(0))    exit q  end if  i = queue(queue_head)  queue_head = mod(queue_head, size(graph, 1) + 1) + 1  do j = 1, size(graph, 1)    if (.not. visited(j)) then      if (graph(j, i) - flow(j, i) > eps) then        visited(j) = .true.        prev(j) = i        if (j == 2) then          <...>          exit q        end if        queue(queue_tail) = j        queue_tail = mod(queue_tail, size(graph, 1) + 1) + 1      end if    end if  end doend do q

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

То, что между узлами и есть ребро, через которое можно пропустить поток, записывается как graph(j, i) - flow(j, i) > eps.

Осталось собрать путь, используя информацию из prev:

if (j == 2) then  temp_capacity = 8  allocate(path(temp_capacity))  temp_len = 0  k = j  do    temp_len = temp_len + 1    if (temp_len > temp_capacity) then      temp_capacity = 2 * temp_capacity      allocate(temp(temp_capacity))      temp(1:temp_len - 1) = path(1:temp_len - 1)      call move_alloc(temp, path)    end if    path(temp_len) = k    k = prev(k)    if (k == 0) exit  end do  if ( temp_len < temp_capacity ) then    allocate(temp(temp_len))    temp(1:temp_len) = path(1:temp_len)    call move_alloc(temp, path)  end if  exit qend if

Здесь используется встроенная процедура move_alloc, которая передаёт владение выделенной памятью от первого аргумента второму.

Полностью поиск кратчайшего пути выглядит так:

function find_shortest_path(graph, flow) result(path)  use, intrinsic :: iso_fortran_env, only: dp => real64  real(dp), intent(in), dimension(:, :) :: graph  real(dp), intent(in), dimension(:, :) :: flow  integer, dimension(:), allocatable :: path  real(dp), parameter :: eps = 1.0e-10_dp  integer, dimension(size(graph, 1) + 1) :: queue  integer :: queue_head, queue_tail  logical, dimension(size(graph, 1)) :: visited  integer, dimension(size(graph, 1)) :: prev  integer :: i, j  integer :: temp_capacity  integer :: temp_len  integer, dimension(:), allocatable :: temp  integer :: k  if (size(graph, 1) /= size(graph, 2)) error stop "hidden_flow:find_shortest_path:arg"  if (size(flow, 1) /= size(flow, 2)) error stop "hidden_flow:find_shortest_path:arg"  if (size(flow, 1) /= size(graph, 1)) error stop "hidden_flow:find_shortest_path:arg"  if (size(graph, 1) < 2) error stop "hidden_flow:find_shortest_path:arg"  queue(1) = 1  queue_head = 1  queue_tail = 2  visited = .false.  visited(1) = .true.  prev(1) = 0  q: do    if (queue_tail == queue_head) then      allocate(path(0))      exit q    end if    i = queue(queue_head)    queue_head = mod(queue_head, size(graph, 1) + 1) + 1    do j = 1, size(graph, 1)      if (.not. visited(j)) then        if (graph(j, i) - flow(j, i) > eps) then          visited(j) = .true.          prev(j) = i          if (j == 2) then            temp_capacity = 8            allocate(path(temp_capacity))            temp_len = 0            k = j            do              temp_len = temp_len + 1              if (temp_len > temp_capacity) then                temp_capacity = 2 * temp_capacity                allocate(temp(temp_capacity))                temp(1:temp_len - 1) = path(1:temp_len - 1)                call move_alloc(temp, path)              end if              path(temp_len) = k              k = prev(k)              if (k == 0) exit            end do            if ( temp_len < temp_capacity ) then              allocate(temp(temp_len))              temp(1:temp_len) = path(1:temp_len)              call move_alloc(temp, path)            end if            exit q          end if          queue(queue_tail) = j          queue_tail = mod(queue_tail, size(graph, 1) + 1) + 1        end if      end if    end do  end do qend function

Теперь нетрудно написать определение максимального потока методом Эдмондса — Карпа:

function edmonds_karp(graph) result(flow)  use, intrinsic :: iso_fortran_env, only: dp => real64  real(dp), intent(in), dimension(:, :) :: graph  real(dp), dimension(size(graph, 1), size(graph, 1)) :: flow  integer, dimension(:), allocatable :: path  integer :: i  real(dp) :: f  if (size(graph, 1) /= size(graph, 2)) error stop "hidden_flow:edmonds_karp:arg"  if (size(graph, 1) < 2) error stop "hidden_flow:edmonds_karp:arg"  flow = 0.0_dp  do    path = find_shortest_path(graph, flow)    if (size(path) < 2) exit    f = graph(path(size(path) - 1), path(size(path))) - flow(path(size(path) - 1), path(size(path)))    do i = size(path) - 2, 1, -1      f = min(f, graph(path(i), path(i + 1)) - flow(path(i), path(i + 1)))    end do    do i = size(path) - 1, 1, -1      flow(path(i), path(i + 1)) = flow(path(i), path(i + 1)) + f      flow(path(i + 1), path(i)) = flow(path(i + 1), path(i)) - f    end do  end doend function

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

Теперь можно вычислить матрицы перехода:

function calc_step_matrix(psi, u) result(m)  use, intrinsic :: iso_fortran_env, only: dp => real64  complex(dp), intent(inout), dimension(:) :: psi  complex(dp), intent(in), dimension(:, :) :: u  real(dp), dimension(size(psi), size(psi)) :: m  real(dp), dimension(size(psi)) :: start  real(dp), dimension(size(psi)) :: finish  real(dp), dimension(size(psi), size(psi)) :: through  real(dp), dimension(size(psi) * 2 + 2, size(psi) * 2 + 2) :: graph  real(dp), dimension(size(psi) * 2 + 2, size(psi) * 2 + 2) :: flow  integer :: i, j  real(dp), parameter :: eps = 1.0e-10_dp  if (size(psi) /= size(u, 1)) error stop "hidden_flow:calc_step_matrix:arg"  if (size(psi) /= size(u, 2)) error stop "hidden_flow:calc_step_matrix:arg"  start = abs(psi) ** 2  psi = matmul(u, psi)  finish = abs(psi) ** 2  through = abs(u)  graph = build_flow_graph(start, finish, through)  flow = edmonds_karp(graph)  do i = 1, size(psi)   do j = 1, size(psi)     if (start(i) <= eps) then       m(j, i) = 0.0_dp     else       m(j, i) = flow(2 + size(psi) + j, 2 + i) / start(i)     end if   end do  end doend function

Отдельно обработывается случай, если $|\psi^{(t-1)}_i|^2 = 0$ : на ноль лучше лишний раз не делить. В таком случае частица не может оказаться в соответствующем состоянии, так что не столь важно как заполнить соответствующую колонку (она же $\Gamma^{(t)}$ ), поэтому мы просто заполняем её нулями.

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

function calc_steps(psi, u, conf) result(path)  use, intrinsic :: iso_fortran_env, only: dp => real64, i8 => int8  complex(dp), intent(inout), dimension(:) :: psi  complex(dp), intent(in), dimension(:, :, :) :: u  integer(i8), intent(in) :: conf  integer(i8), dimension(size(u, 3) + 1) :: path  real(dp), dimension(size(psi), size(psi)) :: m  integer :: i  integer(i8) :: j  real(dp) :: p  real(dp) :: q  if (size(u, 1) /= size(psi)) error stop "hidden_flow:calc_steps:arg"  if (size(u, 2) /= size(psi)) error stop "hidden_flow:calc_steps:arg"      path(1) = conf  do i = 1, size(u, 3)    m = calc_step_matrix(psi, u(:, :, i))    call random_number(p)    path(i + 1) = size(psi) ! prevent round error strike    q = 0.0_dp    do j = 1, size(psi)      q = q + m(j, path(i))      if (q > p) then        path(i + 1) = j        exit      end if    end do  end doend function

Ну и давайте наконец посчитаем траекторию для нашего мысленного эксперимента. У нас будет всего два состояния (n = 2) и три момента времени (t = 3), то есть два шага. На первом шаге работает первое полупрозрачное зеркало, на втором второе. Зеркала одинаковые, поэтому u_1 и u_2 в коде заданы одинаковыми.

subroutine calc_path  use, intrinsic :: iso_fortran_env, only: dp => real64, i8 => int8  integer, parameter :: n = 2  complex(dp), dimension(n), parameter :: psi_0 = [(1.0_dp, 0.0_dp), (0.0_dp, 0.0_dp)]  real(dp), parameter :: one_over_sqrt_two = 1.0_dp / sqrt(2.0_dp)  complex(dp), dimension(n, n), parameter :: u_1 = reshape([ &    -one_over_sqrt_two, one_over_sqrt_two, &    one_over_sqrt_two, one_over_sqrt_two &  ], [n, n])  complex(dp), dimension(n, n), parameter :: u_2 = reshape([ &    -one_over_sqrt_two, one_over_sqrt_two, &    one_over_sqrt_two, one_over_sqrt_two &  ], [n, n])  integer, parameter :: t = 3  integer(i8), dimension(t) :: path  complex(dp), dimension(n) :: psi  complex(dp), dimension(n, n, t - 1) :: u  u(:, :, 1) = u_1  u(:, :, 2) = u_2  psi = psi_0  path = calc_steps(psi, u, 1_i8)  print *, pathend subroutine

Запустим программу несколько раз и увидим то такой вывод:

    1    2    1

то такой:

    1    1    1

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

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

ссылка на оригинал статьи https://habr.com/ru/articles/1099082/