хочу научиться компьютерным экспериментам в математике
текущий план — попробовать записывать сюда происходящее (прямо примеры готового кода и т.п.) не реже раза в неделю
начиная с чего-то совсем простого, а там посмотрим, как пойдет
--Г.М.
текущий план — попробовать записывать сюда происходящее (прямо примеры готового кода и т.п.) не реже раза в неделю
начиная с чего-то совсем простого, а там посмотрим, как пойдет
--Г.М.
def prime_factors(num):
factor = 2
while num > 1:
if num % factor == 0:
yield factor
num = num // factor
else:
factor += 1
for k in [2,3,5,7,11,13,17,19,23,29]:
print("2^",k,"-1",end="=",sep='')
print(*prime_factors(2**k-1),sep='x')
в статью в Квантике про числа Мерсенна (вошла в №1 за 2025 год) захотелось включить короткую программу, которая демонстрирует, что числа вида 2^p-1 не всегда являются простыми
2^2-1=3
2^3-1=7
2^5-1=31
2^7-1=127
2^11-1=23x89
2^13-1=8191
2^17-1=131071
2^19-1=524287
2^23-1=47x178481
2^29-1=233x1103x2089
e^x = 1+x+x²/2+x³/3!+…
пусть x какое-нибудь конкретное число (на картинке x=20) — какие слагаемые дают основной вклад в сумму?
нетрудно понять, что для e^n самое большое слагаемое — это n^n/n!
дальше слагаемые уменьшаются как минимум в геометрической прогрессии и можно оценить
n^n/n! < e^n < (2n+1) n^n/n!
что дает простое доказательство слабой формы формулы Стирлинга:
(n/e)^n < n! < Cn (n/e)^n
прикольно, что это не так далеко от правильной асимптотики n! ~ C√n (n/e)^n
наверное оценку можно усилить, если понять какой ширины «горб» на картинке (в духе цпт и т.п.), но можно ли сделать это не опираясь на Стирлинга? не продумывал еще
===
нравится, что можно ставить компьютерные эксперименты без всякого программирования — в Экселе получается и доступно, и наглядно
файл в комментариях
пусть x какое-нибудь конкретное число (на картинке x=20) — какие слагаемые дают основной вклад в сумму?
нетрудно понять, что для e^n самое большое слагаемое — это n^n/n!
дальше слагаемые уменьшаются как минимум в геометрической прогрессии и можно оценить
n^n/n! < e^n < (2n+1) n^n/n!
что дает простое доказательство слабой формы формулы Стирлинга:
(n/e)^n < n! < Cn (n/e)^n
прикольно, что это не так далеко от правильной асимптотики n! ~ C√n (n/e)^n
наверное оценку можно усилить, если понять какой ширины «горб» на картинке (в духе цпт и т.п.), но можно ли сделать это не опираясь на Стирлинга? не продумывал еще
===
нравится, что можно ставить компьютерные эксперименты без всякого программирования — в Экселе получается и доступно, и наглядно
файл в комментариях
Не планировал еще одного поста так скоро, но пришел Костя Кноп и рассказал забавное:
Как быстро посчитать много знаков числа π, если есть калькулятор с тригонометрическими функциями (но без обратных тригонометрических)? Предлагается такой рецепт: начинаем с грубого приближения π≈3, а дальше делаем пару-тройку итераций x→x+sin(x).
Можно было бы не писать никакой программы — но раз уж тут канал про компьютерные эксперименты:
Кто запустил код — может видеть, что через 4 итерации ошибка уже порядка 10^{-100} (!)
Ну и объяснить это не сложно: за итерацию мы заменяем π+t наπ+t-sin(t), а t-sin(t)≈t³/6 — вот и увеличивается за каждую итерацию количество правильных цифр примерно втрое .
Еще можно заметить, что мы недалеко ушли от методаНьютона (буквально метод Н. для sin(x) — это x→x-tg(x), но рядом с π это примерно то же) .
Тему вычисления знаков π и т.п. наверное еще продолжим.
===
Я ничего не понимаю в библиотеках питона, а mpmath — то что сходу нашлось в гугле по моему запросу (чтобы можно было много цифр синуса считать и т.п.).
Как быстро посчитать много знаков числа π, если есть калькулятор с тригонометрическими функциями (но без обратных тригонометрических)? Предлагается такой рецепт: начинаем с грубого приближения π≈3, а дальше делаем пару-тройку итераций x→x+sin(x).
Можно было бы не писать никакой программы — но раз уж тут канал про компьютерные эксперименты:
from mpmath import *
mp.dps = 200
def iterate(approx):
return approx+sin(approx)
mypi = 3
for k in range(4):
mypi = iterate(mypi)
diff = pi-mypi
print(k+1,":",nstr(mypi,50),
"diff:",nstr(diff,2))
Кто запустил код — может видеть, что через 4 итерации ошибка уже порядка 10^{-100} (!)
Ну и объяснить это не сложно: за итерацию мы заменяем π+t на
Еще можно заметить, что мы недалеко ушли от метода
Тему вычисления знаков π и т.п. наверное еще продолжим.
===
Я ничего не понимаю в библиотеках питона, а mpmath — то что сходу нашлось в гугле по моему запросу (чтобы можно было много цифр синуса считать и т.п.).
набросал вчера относительно длинный черновик поста с примерами того, что хотел бы закодить, но не могу (и почему именно не могу…)
но сегодня что-то посмотрел на один из примеров и решил, что глаза боятся, а… короче, код получился недлинный, но унесу все ж в комментарии
===
интересно считать количество разбиений разных фигур на доминошки
про прямоугольники 2×N и числа Фибоначчи и так все знают
а вот на картинке ацтекский бриллиант порядка 4 — сколькими способами его можно разбить на доминошки? а бриллиант порядка N?
тут общий ответ можно угадать по ответам для небольших N, но перебирать руками тяжеловато (для N=4 вариантов уже больше тысячи), компьютерный эксперимент действительно полезен
про разные способы этот ответ доказывать есть интересная брошюра Е.Смирнова «Три взгляда на ацтекский бриллиант», https://mccme.ru/free-books/dubna/smirnov-aztec.pdf
***
более естественная, казалось бы, идея — считать замощения обычного квадрата N×N… но там по ответам для маленьких N ничего угадать не возможно
все же перебрал на компутере и замощения квадрата 8×8 — в качестве ответа получилось название лекции С.Смирнова https://www.mathnet.ru/present18250
но сегодня что-то посмотрел на один из примеров и решил, что глаза боятся, а… короче, код получился недлинный, но унесу все ж в комментарии
===
интересно считать количество разбиений разных фигур на доминошки
про прямоугольники 2×N и числа Фибоначчи и так все знают
а вот на картинке ацтекский бриллиант порядка 4 — сколькими способами его можно разбить на доминошки? а бриллиант порядка N?
тут общий ответ можно угадать по ответам для небольших N, но перебирать руками тяжеловато (для N=4 вариантов уже больше тысячи), компьютерный эксперимент действительно полезен
про разные способы этот ответ доказывать есть интересная брошюра Е.Смирнова «Три взгляда на ацтекский бриллиант», https://mccme.ru/free-books/dubna/smirnov-aztec.pdf
***
более естественная, казалось бы, идея — считать замощения обычного квадрата N×N… но там по ответам для маленьких N ничего угадать не возможно
все же перебрал на компутере и замощения квадрата 8×8 — в качестве ответа получилось название лекции С.Смирнова https://www.mathnet.ru/present18250
как получить случайное разбиение фигуры на доминошки?
можно попробовать начать с какого-то фиксированного разбиения (например, со всеми горизонтальными доминошками), а дальше применять случайные флипы: выбираем случайный квадрат 2×2 в фигуре, если он разбит на две доминошки, делаем флип (= <-> ||)
(задача: 1) доказать, что флипами можно дойти от любого разбиения квадрата на доминошки до любого другого; 2) придумать фигуру, для разбинений которой на доминошки такое неверно)
на первой картинке для квадрата 20×20 без угла сделано 100 флипов и видно, что мы пока от начального разбиения недалеко ушли; на второй картинке 10000 флиплов — и это выглядит уже довольно хаотично
вроде все ожидаемо, ничего особо интересного
...
попробуем теперь то же самое сделать с ацтекским брильянтом порядка 20… ну первый сюрприз, что после 10 тыс. флипов не происходит почти ничего… это, впрочем, не удивительно, в начале даже подходящих квадратов совсем мало…
картинки с брильянтами и продолжение — в следующем посте
===
способ хранения разбиения из предыдущего поста оказался неудобным тем, что значок не говорит, в какую сторону продолжается покрывающая доминошку клетка — стал хранить в каждой клетке указатель на парную к ней (не в программистском смысле указатель, а один из символов →←↓↑ )
картинки рисую matplotlib'ом (код в комментариях)
каждая доминошка покрывает ровно одну клетку с четной суммой координат — цвет доминошки соответствует тому, в какую сторону из этой клетки продолжается доминошка
можно попробовать начать с какого-то фиксированного разбиения (например, со всеми горизонтальными доминошками), а дальше применять случайные флипы: выбираем случайный квадрат 2×2 в фигуре, если он разбит на две доминошки, делаем флип (= <-> ||)
(задача: 1) доказать, что флипами можно дойти от любого разбиения квадрата на доминошки до любого другого; 2) придумать фигуру, для разбинений которой на доминошки такое неверно)
на первой картинке для квадрата 20×20 без угла сделано 100 флипов и видно, что мы пока от начального разбиения недалеко ушли; на второй картинке 10000 флиплов — и это выглядит уже довольно хаотично
вроде все ожидаемо, ничего особо интересного
...
попробуем теперь то же самое сделать с ацтекским брильянтом порядка 20… ну первый сюрприз, что после 10 тыс. флипов не происходит почти ничего… это, впрочем, не удивительно, в начале даже подходящих квадратов совсем мало…
картинки с брильянтами и продолжение — в следующем посте
===
способ хранения разбиения из предыдущего поста оказался неудобным тем, что значок не говорит, в какую сторону продолжается покрывающая доминошку клетка — стал хранить в каждой клетке указатель на парную к ней (
картинки рисую matplotlib'ом (код в комментариях)
каждая доминошка покрывает ровно одну клетку с четной суммой координат — цвет доминошки соответствует тому, в какую сторону из этой клетки продолжается доминошка
на картинке — генерация случайного разбиения ацтекского брильянта порядка 40 на доминошки (10 тыс. попыток флипов, 100 тыс., 1 млн., много млн¹)
начал с «большого числа» 10000 итераций, а когда дошел до 1000000 — стало казаться, что ничего не получится… а на самом деле это уже почти успех, важно только не опускать руки
может показаться, что и последняя картинка еще недостаточно перемешана — что это за одноцветные области в углах?! но дальше картинка качественно не меняется, бери 10 млн или 100 млн итераций
это проявление замечательной «теоремы о полярном круге» — для случайного разбиения большого брильянта внутри круга хаос, а снаружи всё замерзло в фиксированном состоянии
ну и это еще 40-брильянт не такой уж, на самом деле, большой — в комментарии положу картинку для порядка 100… когда она догенерируется
картинок брильянтов можно найти сколько угодно в интернетах, в «Мат. байках» Витя Клепцын объясняет качественно почему так получается² — но всё равно когда сам запускаю программу и реально этот полярный круг возникает, ощущаю чудо
¹ в середине века у Роллс-Ройсов в техпаспорте в графе «мощность двигателя» было написано просто «достаточная» («adequate»)
² см. https://www.group-telegram.com/mathtabletalks/733 и далее
===
в следующий раз, наверное, будет что-то low tech — мб с экселем вместо питона
начал с «большого числа» 10000 итераций, а когда дошел до 1000000 — стало казаться, что ничего не получится… а на самом деле это уже почти успех, важно только не опускать руки
может показаться, что и последняя картинка еще недостаточно перемешана — что это за одноцветные области в углах?! но дальше картинка качественно не меняется, бери 10 млн или 100 млн итераций
это проявление замечательной «теоремы о полярном круге» — для случайного разбиения большого брильянта внутри круга хаос, а снаружи всё замерзло в фиксированном состоянии
ну и это еще 40-брильянт не такой уж, на самом деле, большой — в комментарии положу картинку для порядка 100… когда она догенерируется
картинок брильянтов можно найти сколько угодно в интернетах, в «Мат. байках» Витя Клепцын объясняет качественно почему так получается² — но всё равно когда сам запускаю программу и реально этот полярный круг возникает, ощущаю чудо
¹ в середине века у Роллс-Ройсов в техпаспорте в графе «мощность двигателя» было написано просто «достаточная» («adequate»)
² см. https://www.group-telegram.com/mathtabletalks/733 и далее
===
в следующий раз, наверное, будет что-то low tech — мб с экселем вместо питона
Рождественская теорема Ферма говорит про то, какие числа представимы в виде суммы двух квадратов
на картинке — компьютерный эксперимент на эту тему: в экселе легко¹ сгенерировать таблицу для сумм двух квадратов
если N=X²+Y², то 2N=(X+Y)²+(X-Y)², поэтому интересна прежде всего представимость нечетных чисел²
вот под основной таблицей посчитано, какие из нечетных чисел до 20 представимы в виде суммы двух квадратов (встречаются в основной таблице), а какие нет
сразу возникает гипотеза…на самом деле, жизнь устроена несколько сложнее, но истина все же скорее возникает из заблуждения, чем из неясности, и дальше уже есть, куда двигаться
продолжение следует
¹ если нужны подробности, как это сделать — см. комментарии
² вместо алгебры можно думать про геометрию (а вместо экселя можно, конечно, пользоваться питоном) — см., например, статью про косые квадраты в https://kvantik.com/issue/pdf/2021-07.pdf
на картинке — компьютерный эксперимент на эту тему: в экселе легко¹ сгенерировать таблицу для сумм двух квадратов
если N=X²+Y², то 2N=(X+Y)²+(X-Y)², поэтому интересна прежде всего представимость нечетных чисел²
вот под основной таблицей посчитано, какие из нечетных чисел до 20 представимы в виде суммы двух квадратов (встречаются в основной таблице), а какие нет
сразу возникает гипотеза…
продолжение следует
¹ если нужны подробности, как это сделать — см. комментарии
² вместо алгебры можно думать про геометрию (а вместо экселя можно, конечно, пользоваться питоном) — см., например, статью про косые квадраты в https://kvantik.com/issue/pdf/2021-07.pdf