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

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

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

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

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

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

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

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

Исходными данными для расчёта служит начальное квантовое состояние и набор унитарных матриц , которые задают динамику квантового состояния:

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

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

Наша цель — превратить и в матрицы перехода такие, что если на -м шаге вероятность пребывания в состоянии равна , а в состоянии 2 — , то

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

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

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

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

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

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

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 do
end 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) :: queue
integer :: queue_head, queue_tail

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

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

logical, dimension(size(graph, 1)) :: visited
integer, dimension(size(graph, 1)) :: prev

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

integer :: temp_capacity
integer :: temp_len
integer, dimension(:), allocatable :: temp

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

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
          <...>
          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 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 q
end 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 q
end 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 do
end 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 do
end function

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

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

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 do
end 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 *, path
end subroutine

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

    1    2    1

то такой:

    1    1    1

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

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