Какой следующий член…? — Ищем формулу для n-го члена последовательности, производящие функции и Z-преобразование13.11.2019 20:33

vxx1ifegfiwsz3zvakjqvc5mop0.pngСкачать файл с кодом и данные можно в оригинале поста в моем блоге

В языке Wolfram Language есть четыре совершенно потрясающие функции: FindSequenceFunction, RSolve, DifferenceRootReduce и FindFormula. В этой статье мы обсудим их возможности и поговорим о функциях, тесно с ними связанных — для поиска параметров линейной рекурсии FindLinearRecurrence (коэффициентов линейного рекуррентного уравнения), производящих функциях GeneratingFunction и Z-преобразовании ZTransform.

Первая функция — FindSequenceFunction — по последовательности чисел ищет выражение для её n-го члена не требуя вообще ничего более.

Hold @ FindSequenceFunction[{1, 1, 2, 3, 5, 8, 13}, n]


80b87533c88cd0137529ae4f4937c355.png

FindSequenceFunction[
{-2, 4Sqrt[Pi],
-16, 16Sqrt[Pi],
-128/3, 32Sqrt[Pi],
-1024/15, 128Sqrt[Pi]/3,
-8192/105, 128Sqrt[Pi]/3},
n]


f338bed8df81abf3f991d904af894454.png
Вторая функция — RSolve — решает рекуррентные уравнения самых разных типов. Элементы могут иметь вид $a[f[n]]$, $a[f[f[n]]]$, $a[f[f[\text{...}f[n]\text{...}]]]$, где f имеет вид: n+A (арифметические разностные уравнения), B*n — геометрические или q-разностные уравнения), B*n+a (арифметико-геометрические функциональные разностные уравнения), B*n^d (степеные геометрические функциональные разностные уравнения), (A*n+B)/(C*n+D) (линейные дробные функциональные разностные уравнения).

RSolve[
	{
		a[n + 3]==2 * a[n],
		a[1]==α,
		a[2]==β,
		a[3]==γ
	},
	a, n
]


5656b0318b4e86118675ddebcabd5f36.png

RSolve[
	{
		v[n]==(2 * Pi * v[n - 2]) / n,
		v[2]==Pi,
		v[3]==(4 * Pi) / 3
	},
	v @ n, n
]


a6deb391cacf5d8db3f4a90b1c18da55.png

Третья функция — DifferenceRootReduce — ищет рекуррентное соотношение для последовательности чисел, n-й член которой имеет заданный вид.

DifferenceRootReduce[-2 * n * Pi * Factorial[(n * 2) - 1],
	n
]


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
]


fbf1daf36c54e6df3150ebf7c207beed.png

Эта функция может много чего ещё, скажем, проверять тождества относительно последовательностей, к примеру: 

DifferenceRootReduce[Fibonacci[2 * n]==Fibonacci[n] * LucasL[n], n]


369b3ef45a3ad032ada5b0aaf9683d2b.png

Здесь LucasL — последовательность чисел Люка (это, по сути, последовательность Фибоначчи, только первые члены не 1, 1, а 1, 3.

Hold @ DifferenceRootReduce @ LucasL @ n


f8e59bd0afe6b4c5434ddc366cd8173c.png

DifferenceRootReduce[LucasL[n]==Fibonacci[n - 1] + Fibonacci[n + 1]]


369b3ef45a3ad032ada5b0aaf9683d2b.png

Как найти рекуррентную формулу для последовательности?


Метод поиска общего члена последовательности часто основан на том, что нужно подобрать рекуррентное уравнение.

Работать это может примерно так: пусть мы ищем n-й член последовательности в виде $f[n]=\sum_{i=1}^ka[i]f[n-i]$. Пусть у нас есть первые члены последовательности: 

sequence = {1, 0, 1, 2, 5, 12, 29, 70, 169, 408, 985, 2378, 5741, 13860, 33461}


1e0672a23dc37c5ee057f82375785fe0.png

Попробуем найти выражение для n-го члена в виде $f[n]=\sum_{i=1}^1a[i]f[n-i]=a[1]f[n-1]$: 

seauenseEq1 = MovingMap[
	Function[
		Dot[Part[#, 1;;1], {a @ 1}]==Part[#, -1]
	],
	sequence, 1
]


2bc4f58d05e3e7751b103d5354dfa50d.png

Hold @ Solve @ seauenseEq1


48792c896bbad9080fb3a1ff9d2b74d1.png

Как видно, решений нет.

Попробуем искать теперь в виде $f[n]=\sum_{i=1}^2a[i]f[n-i]=a[1]f[n-1]+a[2]f[n-2]$: 

seauenseEq2 = MovingMap[
	Function[
		Dot[Part[#, 1;;2], {a @ 1, a @ 2}]==Part[#, -1]
	],
	sequence, 2
]


f0b6e6e5fef79295fb13328c88524aaf.png

Hold @ Solve @ seauenseEq2


20e0d80f4a15e84ca0657538fb36be79.png

Как видим, получилось. Значит, n-й член имеет вид: $f[n]=f[n-1]+2f[n-2]$.

На самом деле есть встроенная функция FindLinearRecurrence, которая позволяет найти линейную рекурсию, подобно тому, как мы это только что сделали: 

Hold @ FindLinearRecurrence @ sequence


cc57304a7e860ddbd3ed7271084d2aad.png

Используя функцию LinearRecurrence можно продлить последовательность: 

LinearRecurrence[{2, 1}, sequence[[1;;2]], 50]


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]


2def6ebe03dc9d428567f2a2ea4fe76a.png

Hold @ sequenseExtension[{1, 2, 2, 1, 1, 2, 2, 1}, 20]


d505af46c84260833e250919c10a448c.png

Hold @ sequenseExtension[
{1, 0, -1, 0, 2, 0, -2, 0, 3, 0, -3, 0, 4, 0, -4},
25
]


ff26b5a3b37f6098e970dee89a84b731.png

Как найти формулу для n-го члена последовательности?


Z-преобразование


Z-преобразование состоит в вычислении ряда вида $\sum_{n=0}^{\infty}f(n)z^{-n}$ от дискретной функции $f(n)$. Это преобразование позволяет свести рекуррентное уравнение для задания последовательности к уравнению относительно образа функции $f(n)$, что аналогично преобразованию Лапласа, которое сводит дифференциальные уравнения к алгебраическим.

Вот как это работает: 

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
]


5cd7e08ab202f0913b7e7b2bb92e221c.png

Посмотрим на примере, скажем, возьмем хорошо известную последовательность Фибоначчи: 

fibonacciEq = f[n]==f[n - 1] + f[n - 2];

initialConditions = {f[1] -> 1, f[2] -> 1};


Ясно, что её стоит переписать в виде, как показано ниже, чтобы не появлялись конструкции типа $f(-1)$ после применения Z-преобразования.

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]]


12739ae3c0836076ca2f6ce12c267f94.png

Решим уравнение относительно образа функции f — ZTransform[f[n], n, z]: 

fZTransformed = ReplaceAll[
	ZTransform[f @ n, n, z],
	Part[Solve[fibonacciEqZTransformed, ZTransform[f @ n, n, z]], 1]
]


ba0d06e4ee5cac2520d2cbfa9e1c8869.png

Выполним обратное Z-преобразование, подставив одновременно начальные условия (заменим n на n-1 в финальном выражении, чтобы наша последовательность имела правильную индексацию (с первого, а не нулевого члена): 

ReplaceAll[InverseZTransform[fZTransformed /. initialConditions, z, n],
	n -> (n - 1)
]


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
]


d2702ee179ce6c4e17679abfd0034091.png

RSolve[
	{
		f[n + 2]==(2 * f[n + 1]) + -(5 * f[n]),
		f[1]==20,
		f[2]==0
	},
	f, n
]


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
]


f1e8a89e313792d886adc87b67921121.png

a79da4ae9770f90b8a29ab135eaa84f5.png

7b7526917cbe1e8f375dd118a22a594a.png

Производящие функции


Производящая функция последовательности $a(n)$ это такая функция $G(x)$, разложение которой в ряд Тейлора (или, более широко, Лорана) имеет вид — $G(x)=\sum_{i=0}^{\infty}a(n)x^n$. Другими словами, коэффициенты при степенях x в разложении функции в ряд задают нашу последовательность.

Скажем, функция $G(x)=\frac{1}{1-x}$ является производящей функцией последовательности 1, 1, 1, 1, …: 

Series[1 / (1 + -x), {x, 0, 10}]


abf7af9eb15f145073a033f5f2ee912d.png

А функция $G(x)=\frac{1}{1-x-x^2}$ является производящей функцией последовательности Фибоначчи 1, 1, 2, 3, 5, 8, 13, …: 

Series[(1 * 1) + (-x) + -(x * 2),
	{x, 0, 10}
]


218fa8e132fb027d073adff2c68f4165.png

Ещё есть разновидность производящей функции — экспоненциальная производящая функция, которая для последовательности $a(n)$ имеет вид — $G(x)=\sum_{i=0}^{\infty}\frac{a(n)}{n!}x^n$.

Скажем, для последовательностей 1, 1, 1, 1… и 1, 1, 2, 3, 5, 8, 13, … экспоненциальные производящие функции таковы — $e^x$ и $\frac{1}{\sqrt{5}}e^{-\frac{2x}{1+\sqrt{5}}}\left(e^{\sqrt{5}x}-1\right)$: 

ReplaceAll[Normal[Series[E ^ x, {x, 0, 10}]],
	Power[x, n_] :> ((x ^ n) * Factorial[n])
]


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])
]


e18e27905aab6bc5d1a090fb6ae305f9.png

Производящую функцию в Wolfram Language можно найти двумя функциями — GeneratingFunction и FindGeneratingFunction (экспоненциальную с помощью ExponentialGeneratingFunction): 

GeneratingFunction[-(m * Factorial[n]), {n, m}, {x, y}]


26cd71180a866543ffd730c1e7368f53.png

TraditionalForm[
	FullSimplify[
		ExponentialGeneratingFunction[-(n * Factorial[n - 1] * Factorial[2 * n]), n, x]
	]
]


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]
]


19c79ccd9b984b853e31f85f277e82b5.png

d7561a061aa1a33e437397ee78658972.png

bf6813c2a7cc3c36e3f0eea54471263e.png

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"


911876cd2344e2f74361e048bdee7f20.png

ResourceFunction[«OEISSequenceData»] — выдает датасет с полной информацией из базы: 

sequenceData[666] = OEISSequenceData[666, "Dataset"]


1757e3c9b1f21aa10f67114c5987fa83.png

Скажем, можно «вытащить» код на языке Wolfram Language: 

Hold @ Normal @ sequenceData[666]["CodeWolframLanguageStrings"]


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}}
	]
]


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]


8e278425af83dc452b3220a45a8b7db2.png

formulas = FindFormula[data, x]


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"
]


615e2819fedf142d8c3eb73ead640f91.png

Можно построить и большее количество зависимостей, скажем, 10: 

formulas = FindFormula[data, x, 10]


d9669df48c3721bac0a32ac0db81c0b1.png

Plot[formulas,
	{x, -10, 10},
	PlotStyle -> AbsoluteThickness[3],
	Prolog -> {AbsolutePointSize[5], LightGray, Point @ data},
	Background -> White, ImageSize -> 800, PlotLegends -> "Expressions"
]


826b4cdc6fbae9753578434b3d719854.png

Стоит отметить, что есть функция, аналогичная по функционалу, которая ищет вероятностное распределение — FindDistribution.

Для сотрудничества — пишите личное сообщение на Хабре или в мою группу ВКонтакте.
Канал YouTube — вебинары и обучающие ролики.
Регистрация на новые курсы. Готовый онлайн курс.

© Habrahabr.ru