![](https://habrastorage.org/webt/vx/x1/if/vxx1ifegfiwsz3zvakjqvc5mop0.png)
В языке Wolfram Language есть четыре совершенно потрясающие функции:
FindSequenceFunction
, RSolve
, DifferenceRootReduce
и FindFormula
. В этой статье мы обсудим их возможности и поговорим о функциях, тесно с ними связанных — для поиска параметров линейной рекурсии FindLinearRecurrence
(коэффициентов линейного рекуррентного уравнения), производящих функциях GeneratingFunction
и Z-преобразовании ZTransform
.Первая функция — FindSequenceFunction — по последовательности чисел ищет выражение для её n-го члена не требуя вообще ничего более.
Hold @ FindSequenceFunction[{1, 1, 2, 3, 5, 8, 13}, n]
![](https://habrastorage.org/getpro/habr/post_images/80b/875/33c/80b87533c88cd0137529ae4f4937c355.png)
FindSequenceFunction[
{-2, 4Sqrt[Pi],
-16, 16Sqrt[Pi],
-128/3, 32Sqrt[Pi],
-1024/15, 128Sqrt[Pi]/3,
-8192/105, 128Sqrt[Pi]/3},
n]
![](https://habrastorage.org/getpro/habr/post_images/f33/8be/d8d/f338bed8df81abf3f991d904af894454.png)
Вторая функция — RSolve — решает рекуррентные уравнения самых разных типов. Элементы могут иметь вид
RSolve[
{
a[n + 3]==2 * a[n],
a[1]==α,
a[2]==β,
a[3]==γ
},
a, n
]
![](https://habrastorage.org/getpro/habr/post_images/565/6b0/318/5656b0318b4e86118675ddebcabd5f36.png)
RSolve[
{
v[n]==(2 * Pi * v[n - 2]) / n,
v[2]==Pi,
v[3]==(4 * Pi) / 3
},
v @ n, n
]
![](https://habrastorage.org/getpro/habr/post_images/a6d/eb3/91c/a6deb391cacf5d8db3f4a90b1c18da55.png)
Третья функция — DifferenceRootReduce — ищет рекуррентное соотношение для последовательности чисел, n-й член которой имеет заданный вид.
DifferenceRootReduce[-2 * n * Pi * Factorial[(n * 2) - 1],
n
]
![](https://habrastorage.org/getpro/habr/post_images/694/dd8/197/694dd81973fd06c87eb3ce5c0b013c4f.png)
RSolve[
{
(-8 * y[n]) + n * y[2 + n]==0,
y[-1]==1/4,
y[0]==0,
y[1]==-2,
y[2]==4Sqrt[Pi]
},
y, n
]
![](https://habrastorage.org/getpro/habr/post_images/fbf/1da/f36/fbf1daf36c54e6df3150ebf7c207beed.png)
Эта функция может много чего ещё, скажем, проверять тождества относительно последовательностей, к примеру:
DifferenceRootReduce[Fibonacci[2 * n]==Fibonacci[n] * LucasL[n], n]
![](https://habrastorage.org/getpro/habr/post_images/369/b3e/f45/369b3ef45a3ad032ada5b0aaf9683d2b.png)
Здесь LucasL — последовательность чисел Люка (это, по сути, последовательность Фибоначчи, только первые члены не 1, 1, а 1, 3.
Hold @ DifferenceRootReduce @ LucasL @ n
![](https://habrastorage.org/getpro/habr/post_images/f8e/59b/d0a/f8e59bd0afe6b4c5434ddc366cd8173c.png)
DifferenceRootReduce[LucasL[n]==Fibonacci[n - 1] + Fibonacci[n + 1]]
![](https://habrastorage.org/getpro/habr/post_images/369/b3e/f45/369b3ef45a3ad032ada5b0aaf9683d2b.png)
Как найти рекуррентную формулу для последовательности?
Метод поиска общего члена последовательности часто основан на том, что нужно подобрать рекуррентное уравнение.
Работать это может примерно так: пусть мы ищем n-й член последовательности в виде
sequence = {1, 0, 1, 2, 5, 12, 29, 70, 169, 408, 985, 2378, 5741, 13860, 33461}
![](https://habrastorage.org/getpro/habr/post_images/1e0/672/a23/1e0672a23dc37c5ee057f82375785fe0.png)
Попробуем найти выражение для n-го члена в виде
seauenseEq1 = MovingMap[
Function[
Dot[Part[#, 1;;1], {a @ 1}]==Part[#, -1]
],
sequence, 1
]
![](https://habrastorage.org/getpro/habr/post_images/2bc/4f5/8d0/2bc4f58d05e3e7751b103d5354dfa50d.png)
Hold @ Solve @ seauenseEq1
![](https://habrastorage.org/getpro/habr/post_images/487/92c/896/48792c896bbad9080fb3a1ff9d2b74d1.png)
Как видно, решений нет.
Попробуем искать теперь в виде
seauenseEq2 = MovingMap[
Function[
Dot[Part[#, 1;;2], {a @ 1, a @ 2}]==Part[#, -1]
],
sequence, 2
]
![](https://habrastorage.org/getpro/habr/post_images/f0b/6e6/e5f/f0b6e6e5fef79295fb13328c88524aaf.png)
Hold @ Solve @ seauenseEq2
![](https://habrastorage.org/getpro/habr/post_images/20e/0d8/0f4/20e0d80f4a15e84ca0657538fb36be79.png)
Как видим, получилось. Значит, n-й член имеет вид:
На самом деле есть встроенная функция
FindLinearRecurrence
, которая позволяет найти линейную рекурсию, подобно тому, как мы это только что сделали:Hold @ FindLinearRecurrence @ sequence
![](https://habrastorage.org/getpro/habr/post_images/cc5/730/4a7/cc57304a7e860ddbd3ed7271084d2aad.png)
Используя функцию
LinearRecurrence
можно продлить последовательность:LinearRecurrence[{2, 1}, sequence[[1;;2]], 50]
![](https://habrastorage.org/getpro/habr/post_images/d4e/d65/23e/d4ed6523eac9ff4c763545e935e45fa1.png)
Или объединить все в одну строчку, построив функцию, которая: продлит последовательность, выдаст разностное уравнение и найдет общую формулу для n-го члена:
sequenseExtension[list_, n_] := Module[
{lr, eq},
lr = FindLinearRecurrence @ list;
eq = Flatten[
{
a[k]==Total[
Table[
a[k + -i] * Part[lr, i],
{i, 1, Length @ lr}
]
],
Table[a[i], list[[i]]], {i, 1, Length @ lr}]
}
];
<|
"Уравнение" -> eq,
"Формула" -> FullSimplify[a[k] /. Part[RSolve[eq, a, k], 1]],
"Продление" -> LinearRecurrence[lr, Part[list, Span[1, Length[lr]]], n]
|>
];
Hold @ sequenseExtension[{1, 1, 2, 3, 5}, 20]
![](https://habrastorage.org/getpro/habr/post_images/2de/f6e/be0/2def6ebe03dc9d428567f2a2ea4fe76a.png)
Hold @ sequenseExtension[{1, 2, 2, 1, 1, 2, 2, 1}, 20]
![](https://habrastorage.org/getpro/habr/post_images/d50/5af/46c/d505af46c84260833e250919c10a448c.png)
Hold @ sequenseExtension[
{1, 0, -1, 0, 2, 0, -2, 0, 3, 0, -3, 0, 4, 0, -4},
25
]
![](https://habrastorage.org/getpro/habr/post_images/ff2/6b5/a3b/ff26b5a3b37f6098e970dee89a84b731.png)
Как найти формулу для n-го члена последовательности?
Z-преобразование
Z-преобразование состоит в вычислении ряда вида
Вот как это работает:
Grid[
Transpose[
Function[
{
#,
Map[TraditionalForm, Map[FullSimplify, ZTransform[#, n, z]]]
}
][
{
f[n - 2],
f[n - 1],
f @ n,
f[n + 1],
f[n + 2]
}
]
],
Background -> White, Dividers -> All
]
![](https://habrastorage.org/getpro/habr/post_images/5cd/7e0/8ab/5cd7e08ab202f0913b7e7b2bb92e221c.png)
Посмотрим на примере, скажем, возьмем хорошо известную последовательность Фибоначчи:
fibonacciEq = f[n]==f[n - 1] + f[n - 2];
initialConditions = {f[1] -> 1, f[2] -> 1};
Ясно, что её стоит переписать в виде, как показано ниже, чтобы не появлялись конструкции типа
fibonacciEq = f[n + 2]==f[n + 1] + f[n];
initialConditions = {f[0] -> 1, f[1] -> 1};
Осуществим Z-преобразование:
fibonacciEqZTransformed = ReplaceAll[fibonacciEq, pattern:f[__] :> ZTransform[pattern, n, z]]
![](https://habrastorage.org/getpro/habr/post_images/127/39a/e3c/12739ae3c0836076ca2f6ce12c267f94.png)
Решим уравнение относительно образа функции f — ZTransform[f[n],n,z]:
fZTransformed = ReplaceAll[
ZTransform[f @ n, n, z],
Part[Solve[fibonacciEqZTransformed, ZTransform[f @ n, n, z]], 1]
]
![](https://habrastorage.org/getpro/habr/post_images/ba0/d06/e4e/ba0d06e4ee5cac2520d2cbfa9e1c8869.png)
Выполним обратное Z-преобразование, подставив одновременно начальные условия (заменим n на n-1 в финальном выражении, чтобы наша последовательность имела правильную индексацию (с первого, а не нулевого члена):
ReplaceAll[InverseZTransform[fZTransformed /. initialConditions, z, n],
n -> (n - 1)
]
![](https://habrastorage.org/getpro/habr/post_images/265/394/566/265394566f9efec4011542a91933b272.png)
Естестевенно это можно автоматизировать, создав свой аналог RSolve:
myRSolve[eq_, initials_, f_, n_] := Module[
{z, initialsInner, eqZTransformed, fZTransformed},
initialsInner = ReplaceAll[initials, f[x_] :> f[x - 1]];
eqZTransformed = ReplaceAll[eq, pattern:f[__] :> ZTransform[pattern, n, z]];
fZTransformed = ReplaceAll[ZTransform[f @ n, n, z],
Part[Solve[eqZTransformed, ZTransform[f @ n, n, z]], 1]
];
FullSimplify[
InverseZTransform[fZTransformed /. initialsInner, z, n] /. n -> (n - 1)
]
];
myRSolve[
{
f[n + 2]==(2 * f[n + 1]) + -(5 * f[n])
},
{f[1] -> 20, f[2] -> 0},
f, n
]
![](https://habrastorage.org/getpro/habr/post_images/d27/02e/e17/d2702ee179ce6c4e17679abfd0034091.png)
RSolve[
{
f[n + 2]==(2 * f[n + 1]) + -(5 * f[n]),
f[1]==20,
f[2]==0
},
f, n
]
![](https://habrastorage.org/getpro/habr/post_images/006/6b4/491/0066b44910be89c2964f2d901ad4a02d.png)
Но, конечно, RSolve содержит намного больше возможностей для решения самых разных дискретных уравнений, на которых мы не будем останавливаться подробнее:
RSolve[a[n]==(n * a[n]) + n, a, n],
RSolve[
{
a[n + 1]==(2 * a[n]) + (3 * a[n]) + 4,
a[0]==0
},
a, n
],
RSolve[
y[n + 1 * 3]==(2 * y[n + 1 * 6]) + n * 2,
y, n
]
![](https://habrastorage.org/getpro/habr/post_images/f1e/8a8/9e3/f1e8a89e313792d886adc87b67921121.png)
![](https://habrastorage.org/getpro/habr/post_images/a79/da4/ae9/a79da4ae9770f90b8a29ab135eaa84f5.png)
![](https://habrastorage.org/getpro/habr/post_images/7b7/526/917/7b7526917cbe1e8f375dd118a22a594a.png)
Производящие функции
Производящая функция последовательности
Скажем, функция
Series[1 / (1 + -x), {x, 0, 10}]
![](https://habrastorage.org/getpro/habr/post_images/abf/7af/9eb/abf7af9eb15f145073a033f5f2ee912d.png)
А функция
Series[(1 * 1) + (-x) + -(x * 2),
{x, 0, 10}
]
![](https://habrastorage.org/getpro/habr/post_images/218/fa8/e13/218fa8e132fb027d073adff2c68f4165.png)
Ещё есть разновидность производящей функции — экспоненциальная производящая функция, которая для последовательности
Скажем, для последовательностей 1, 1, 1, 1… и 1, 1, 2, 3, 5, 8, 13,… экспоненциальные производящие функции таковы —
ReplaceAll[Normal[Series[E ^ x, {x, 0, 10}]],
Power[x, n_] :> ((x ^ n) * Factorial[n])
]
![](https://habrastorage.org/getpro/habr/post_images/90a/4ff/4a2/90a4ff4a2c36601893893a8419b52c10.png)
ReplaceAll[
Normal[
FullSimplify[
Series[
Plus[E,
(-(2 * x * 1)) + 5 * ((E * 5 * x) - 1) * 5
],
{x, 0, 10}
]
]
],
Power[x, n_] :> ((x ^ n) * Factorial[n])
]
![](https://habrastorage.org/getpro/habr/post_images/e18/e27/905/e18e27905aab6bc5d1a090fb6ae305f9.png)
Производящую функцию в Wolfram Language можно найти двумя функциями —
GeneratingFunction
и FindGeneratingFunction
(экспоненциальную с помощью ExponentialGeneratingFunction
):GeneratingFunction[-(m * Factorial[n]), {n, m}, {x, y}]
![](https://habrastorage.org/getpro/habr/post_images/26c/d71/180/26cd71180a866543ffd730c1e7368f53.png)
TraditionalForm[
FullSimplify[
ExponentialGeneratingFunction[-(n * Factorial[n - 1] * Factorial[2 * n]), n, x]
]
]
![](https://habrastorage.org/getpro/habr/post_images/dc6/14f/d18/dc614fd182386773d3fb84d4242cafa6.png)
Есть много методов поиска общего члена последовательности с помощью производящих функций. Не будем подробно останавливаться на этом, скажем, только что неплохая теория есть на сайте genfunc.ru.
Один из методов похож на Z-преобразование:
generatingFEq = ReplaceAll[
f[n + 2]==f[n + 1] + f[n],
pattern:f[__] :> GeneratingFunction[pattern, n, z]
],
generatingF = ReplaceAll[
GeneratingFunction[f @ n, n, z],
Part[Solve[generatingFEq, GeneratingFunction[f @ n, n, z]], 1]
],
nthTerm = SeriesCoefficient[generatingF, {z, 0, n}],
FullSimplify[
ReplaceAll[ReplaceAll[nthTerm, {f[0] -> 1, f[1] -> 1}],
n -> (n - 1)
],
GreaterEqual[n, 1]
]
![](https://habrastorage.org/getpro/habr/post_images/19c/79c/cd9/19c79ccd9b984b853e31f85f277e82b5.png)
![](https://habrastorage.org/getpro/habr/post_images/d75/61a/061/d7561a061aa1a33e437397ee78658972.png)
![](https://habrastorage.org/getpro/habr/post_images/bf6/813/c2a/bf6813c2a7cc3c36e3f0eea54471263e.png)
![](https://habrastorage.org/getpro/habr/post_images/265/394/566/265394566f9efec4011542a91933b272.png)
OEIS — Онлайн-энциклопедия целочисленных последовательностей и интеграция с Wolfram Language
В интернете доступна совершенно потрясающая коллекция числовых последовательностей — OEIS (On-Line Encyclopedia of Integer Sequences). Она была создана Нилом Слоуном во время его исследовательской деятельности в AT&T Labs. В OEIS хранится информация о целочисленных последовательностях, представляющих интерес как для любителей, так и для специалистов в математике, комбинаторике, теории чисел, теории игр, физике, химии, биологии, информатике. На данный момент там собрано 329085 последовательностей. Запись в OEIS включает в себя первые элементы последовательности, ключевые слова, математическое описание, фамилии авторов, ссылки на литературу; присутствует возможность построения графика или проигрывания музыкального представления последовательности. Поиск в базе данных может осуществляться по ключевым словам и по подпоследовательности.
Недавно появилась интеграция с этой базой внутри Wolfram Language (при использовании важно понимать, что это разработка пользователей — с недавного времени можно выгружать свой код в репозиторий Wolfram Function Repository). Достаточно просто указать номер интересующей вас последовательности или список номеров.
OEISSequenceData = ResourceFunction @ "OEISSequenceData";
OEISSequence = ResourceFunction @ "OEISSequence";
ResourceFunction[«OEISSequence»] — просто выдает первые члены последовательности:
Hold @ OEISSequence @ "A666"
![](https://habrastorage.org/getpro/habr/post_images/911/876/cd2/911876cd2344e2f74361e048bdee7f20.png)
ResourceFunction[«OEISSequenceData»] — выдает датасет с полной информацией из базы:
sequenceData[666] = OEISSequenceData[666, "Dataset"]
![](https://habrastorage.org/getpro/habr/post_images/175/7e3/c9b/1757e3c9b1f21aa10f67114c5987fa83.png)
Скажем, можно «вытащить» код на языке Wolfram Language:
Hold @ Normal @ sequenceData[666]["CodeWolframLanguageStrings"]
![](https://habrastorage.org/getpro/habr/post_images/175/7e3/c9b/1757e3c9b1f21aa10f67114c5987fa83.png)
Или набор случайно выбранных последовательностей с интересующей по ним информацией:
randomSequences = Dataset @ Map[
Normal,
OEISSequenceData[RandomInteger[{1, 300000}, 10], "Dataset"]
];
Function[
Framed[#, FrameStyle -> None, FrameMargins -> 5, Background -> White]
][
Grid[
Join[
{
Map[Style[#, Bold, 18]&,
{"Название", "Формулы", "Ссылки", "Первые члены", "График первых членов"}
]
},
Map[
Function[
Map[
Function[
TextCell[#, LineIndent -> 0, FontSize -> 12, FontFamily -> "Open Sans Light"]
],
{
Style[Part[#, 1], 16],
Row[Part[#, 4], "\n"],
Row[Part[#, 3], "\n"],
Style[Row[Part[#, 2], "; "], 10],
ListLinePlot[Part[#, 2], ImageSize -> Full]
}
]
],
Values @ Normal @ randomSequences[All, {"Name", "Sequence", "References", "Formulae"}]
]
],
Dividers -> {{None, {LightGray}, None}, {None, {LightGray}, None}},
ItemStyle -> Directive[FontSize -> 12, FontFamily -> "Open Sans Light"],
ItemSize -> {{15, 25, 10, 15, 15}, Automatic},
Alignment -> {Left, Center},
Background -> {None, {LightOrange, White}}
]
]
![](https://habrastorage.org/getpro/habr/post_images/154/4c1/7fd/1544c17fdb6d82599459d556a783d0fe.png)
Поиск потенциально возможной формулы
Наконец, хотелось бы отметить функцию
FindFormula
, которая по заданному набору чисел строит формулу, которая их может описать. Примем зависимостей подобрать можно много и из разных классов функций.data = Table[
{
x,
Sin[2 * x] + Cos[x] + RandomVariate[NormalDistribution[0, 0.2]]
},
{x, RandomReal[{-10, 10}, 1000]}
];
ListPlot[data, Background -> White, ImageSize -> 600]
![](https://habrastorage.org/getpro/habr/post_images/8e2/784/25a/8e278425af83dc452b3220a45a8b7db2.png)
formulas = FindFormula[data, x]
![](https://habrastorage.org/getpro/habr/post_images/1fa/8b1/803/1fa8b1803a27f022add215755ecb750e.png)
Как видно, Wolfram Language подобрал функцию, очень близкую к той, на основе которой были построены «зашумленные» данные, а именно — Sin[2x]+Cos[x]:
Plot[formulas,
{x, -10, 10},
PlotStyle -> AbsoluteThickness[3],
Prolog -> {AbsolutePointSize[5], Gray, Point @ data},
Background -> White, ImageSize -> 800, PlotLegends -> "Expressions"
]
![](https://habrastorage.org/getpro/habr/post_images/615/e28/19f/615e2819fedf142d8c3eb73ead640f91.png)
Можно построить и большее количество зависимостей, скажем, 10:
formulas = FindFormula[data, x, 10]
![](https://habrastorage.org/getpro/habr/post_images/d96/69d/f48/d9669df48c3721bac0a32ac0db81c0b1.png)
Plot[formulas,
{x, -10, 10},
PlotStyle -> AbsoluteThickness[3],
Prolog -> {AbsolutePointSize[5], LightGray, Point @ data},
Background -> White, ImageSize -> 800, PlotLegends -> "Expressions"
]
![](https://habrastorage.org/getpro/habr/post_images/826/b4c/dc6/826b4cdc6fbae9753578434b3d719854.png)
Стоит отметить, что есть функция, аналогичная по функционалу, которая ищет вероятностное распределение —
FindDistribution
.Для сотрудничества — пишите личное сообщение на Хабре или в мою группу ВКонтакте.
Канал YouTube — вебинары и обучающие ролики.
Регистрация на новые курсы. Готовый онлайн курс.