Матрицы и уравнения

Как известно, Льюис Кэрролл не только написал книги про Алису, но и был математиком. Занимался он в основном линейными уравнениями, матрицами и детерминантами. Я о его математических достижениях узнал совсем недавно. Точнее, я прочитал его статью о методе вычисления детерминантов.
Вот что он придумал. Пусть, например, нам нужно вычислить детерминант матрицы

\[\left[\begin{array}{ccccc}4&4&2&1&2\\4&4&5&2&2\\4&2&5&3&1\\2&5&3&1&2\\1&2&0&4&5\end{array}\right]\]

Делаем так: сначала вычисляем детерминанты блоков $2\times 2$ и строим из них матрицу размерности на один меньше:

\[\left[\begin{array}{cccc}0&12&-1&-2\\-8&10&5&-4\\16&-19&-4&5\\-1&-6&12&-3\end{array}\right]\]

Далее делаем то же самое, но теперь каждый детерминант $2\times 2$ делим на соответствующий элемент внутренности первой матрицы. Внутренность матрицы — это матрица, получающаяся из исходной вычеркиванием первой и последней строки и первого и последнего столбца. Получаем

\[\left[\begin{array}{ccc}24&14&7\\-4&11&3\\-23&-84&-48\end{array}\right]\]

Повторяя процесс, получаем

\[\left[\begin{array}{ccc}\begin{array}{cc} 32 & -7 \\-31 & 69\end{array}\end{array}\right]\]

\[\left[181\right]\]

Получившееся число и есть детерминант исходной матрицы.

Если сравнить число операций

\[N(\mbox{Доджсон})=3(n-1)^2+4(n-2)^2+\ldots=\frac{1}{3} (n-1) \left(4 n^2-5 n+3\right)\]

с числом операций при разложении по строке

\[N(\mbox{По строке})=e \Gamma (n+1,1)-2,\]

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

\[N(\mbox{Гаусс})=\frac{1}{6} (n-1) \left(4 n^2+n+6\right)\]

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

" If the determinant of a block $= R$, the determinant of any minor of the $m$th degree of the adjugate block is the product of $R^{m-1}$ and the coefficient which, in $R$, multiplies the determinant of the corresponding minor."

(Jacobi)

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

---

PS Только что выяснил, что описанный метод вычисления детерминантов можно найти в английской википедии в соответствующей статье. Там тождество Доджсона называется Desnanot-Jacobi identity, и это, по-моему, правильно.

Интерференция и SVG

Вчера Ваня напомнил про один вопрос из квантовой механики, с которым я когда-то разобрался.
Все знают, как борновская амплитуда рассеяния на нескольких одинаковых центрах выражается через амплитуду рассеяния на одном. А именно, появляется множитель
\[F(\vec{q})=\sum_{k=0}^{N} e^{i\vec{q}\vec{r}_k}.\]
Если центры образуют периодическую структуру, этот множитель приводит к тому, что в интерференционных максимумах сечение увеличивается в $N^2$ раз.
Допустим теперь, что мы случайно бросаем рассеивающие центры в большой ящик и после этого они остаются неподвижными. Каким будет сечение? Я в свою бытность студентом 3-его курса, пока не обдумал вопрос самостоятельно, считал, что дифференциальное сечение будет примерно равно
\[\frac{d\sigma_N}{d\Omega}\approx N\frac{d\sigma_1}{d\Omega},\]
а отклонения будут расти медленнее чем $N$ (как $\sqrt{N}$). Ниже скрипт, который поможет "экспериментально" изучить вопрос. Как можно заметить, отношение сечения на $N$ центрах к умноженному на $N$ сечению на одном и не думает стремиться к единице при $N\to\infty$. Действительно ли это так, или это артефакт моделирования?

---

$\displaystyle{\frac{d\sigma_N/d\Omega}{N d\sigma_1/d\Omega}=\left|F(q)\right|^2/N=}$?

$\lambda=\frac{2\pi}{q}=$ — расстояние между красными линиями.

---

Вообще-то, этот пост еще служит проверкой возможности включения svg в блог. Скажу честно, делается через ж... использование javascript. Правда, если хочешь интерактивности, javascript все-равно нужен. Кто будет смотреть в код, обратите внимание на функцию createSVG. Вроде, мелочь, а времени много экономит.

Измерения и информация

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

Все, наверное, решали задачу про 12 монет:

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

Задача не тривиальная, но все же, немного подумав, решить можно. Сразу возникает вопрос: а нельзя ли увеличить общее количество монет? Оказывается, что можно до 13. А если есть эталонная, то и до 14. А из скольки монет можно найти фальшивую за четыре взвешивания, за пять?

Сформулируем задачу для любителей обобщений:

Из какого максимального количества монет и каким образом можно найти одну фальшивую за $n$ взвешиваний?

Задачу эту я узнал, на школе Ландау в Черноголовке в далеком 1997 году. На этой школе еще были интересные лекции про криптографию, например, нам рассказали как работает алгоритм шифрования с открытым ключом RSA . Но сейчас не об этом. Со мной был один товарищ, который предложил рекурсивный подход: сначала сгруппировать монеты в кучки, найти среди кучек фальшивую, а затем среди монет этой кучки найти фальшивую. Например, для $n=6$ взвешиваний можно взять аж $13\cdot 14=182$ монеты и сначала применить достижения первой задачи к 13 столбикам по 14 монет, а затем их же к монетам "фальшивого столбика". Хорошая идея, но я решал задачу совсем по-другому и у меня для $6$ взвешиваний получалось $364$ монеты. С тех пор я уже не раз порывался написать какую-нибудь популярную заметку об этой и других подобных задачах, но руки так и не дошли. Я также несколько раз задавал эту задачу коллегам/товарищам, но почему-то у них она не пошла.

Недавно я вспомнил эту историю потому что решал на сайте Diofant.ru другую задачу. Для тех, кому не хочется идти по ссылке, приведу условие:

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

Опять же, в условии можно заменить 7 на $n$. А можно и два радиоактивных (при чем тут радиоактивность?) заменить на $m$.

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

---

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

Скрипт:
Вот еще одна подсказка — рабочий скрипт:

Введите число взвешиваний

 

0		0

 

---

Ответ Максимальное число монет при отсутствии эталонной равно
\[N=\frac{3^n-1}2\]
Максимальное число монет при наличии эталонной равно
\[N=\frac{3^n+1}2\]

Первое взвешивание: $(3^{n-1}-1)/2$ на каждой чашке. Если есть эталонная, то $(3^{n-1}+1)/2$ на одной чашке и $(3^{n-1}-1)/2$ плюс эталонная — на другой. После $k$-ого взвешивания может быть три ситуации:

1. Фальшивая монета находится среди $(3^{n-k}+1)/2$ монет.
2. Фальшивая монета содержится в одной из двух стопок, причем легче она или тяжелее определятся тем, в какой именно стопке она находится. В стопках находится $(3^{n-k}+1)/2$ и $(3^{n-k}-1)/2$ монет (если нужно, можно добавить эталонную).

В зависимости от случая, делаем следующее:

1. В первом случае кладем из стопки под подозрением $(3^{n-k-1}+1)/2$ на одну чашу и $(3^{n-k-1}-1)/2$ плюс эталонная — на вторую. После взвешивания получаем случай 1 или случай 2.
2. Во втором случае на каждую чашу кладем из большей стопки по $(3^{n-k-1}+1)/2$ монет, а из меньшей — по $(3^{n-k-1}-1)/2$. После взвешивания получаем случай 2.

Очевидно, после $k=n$ взвешиваний оба случая однозначно определяют фальшивую.

Ускорение сходимости

Недавно пришлось разобраться в вопросе об ускорении сходимости последовательностей. До этого я был в полной уверенности, что знаю достаточно об этом, а чего не знаю — и не надо никогда. Но оказалось, что бывает надо. Расскажу про один простой метод, описанный в статье Дэвида Бродхерста (Broadhurst). Статья, конечно, не про метод.

Рассмотрим в качестве подопытной свинки $\zeta$-функцию Римана:
\[\zeta(n)=\sum_{k=1}^{\infty}\frac 1{k^n}\]

Просуммировав $N$ членов суммы, получим примерно $P\approx (n-1)\log_{10}N$ десятичных знаков. Можно, конечно, "пришить хвост", т.е., аппроксимировать остаток суммы интегралом, но это не может изменить того факта, что сходимость — логарифмическая $P\propto \log N$. Меняется только коэффициент пропорциональности. Справедливости ради заметим, что для многих практических целей достаточно, чтобы этот коэффициент был достаточно большим. Но что, если нам все-же хочется уметь вычислять такие суммы с высокой точностью?

Замечательный метод ускорения сходимости для таких последовательностей описан в вышеупомянутой статье. Я сначала не особенно сильно на него надеялся. Интуитивно мне казалось, что если мы вычислили $N$ первых членов суммы, то лучшее, что мы с ними можем сделать — просто сложить. Фокус в том, что мы не просто знаем эти члены, но еще и знаем кое-что об остатке $r_k$. А именно, функция $k^{n-1}r_k$ раскладывается в ряд вблизи точки $k=\infty$. Это позволяет построить некоторую взвешенную сумму первых $N$ членов так, что $P\propto N$!

В качестве иллюстрации прилагаю маленький скрипт, написанный по мотивам вышеупомянутой статьи, который вычисляет $\zeta(n)$. Внимательный взгляд найдет пару отличий от описанного в статье. Поскольку написан он на javascript, то точность в какой-то момент теряется (точнее, не растет), но это проблемы вычислений с фиксированной точностью. Вбейте в форму ниже $N$ и $n$ и нажмите кнопку.

---

N= членов суммы для ζ() дают
при тупом суммировании, и
после ускорения.

---

Если кому-то стало интересно, предлагаю выяснить как выглядит вектор весов и определить коэффициент пропорциональности между P и N (экспериментально). Результат может удивить. Еще один вопрос — какой должна быть рабочая точность? Для удобства ниже — упрощенный (без учета потери точности) код скрипта
 function Zeta(N,n) {
  var terms= [0];
  //accumulating from above
  for(var i=N;i>0;i--) {
  terms.push(terms[terms.length-1]+1/Math.pow(i,n));
  };
  //acceleration starts
  for(var m=0;m<N;m++) {
  for(var k=0;k<(N-m);k++) {
  terms[k]=((N-k)*terms[k]-(N-k-m-n+1)*terms[k+1])/(m+n-1)
  };
  };
  //direct sum
  document.getElementById('direct').value=terms[N];
  //accelerated
  document.getElementById('accelerated').value=terms[N]-terms[0];
 }