soft :: текст

Maple: разложение рациональной дроби

Еще один яркий пример использования Maple для алгоритмизации рутинных и часто используемых операций - разложение рациональной дроби с действительными коэффициентами на сумму простейших. Уж более стандартную и столь распространенную задачу как в алгебре, так и в математическом анализе придумать сложно. Существует несколько методов построения сей конструкции, однако остановимся на наиболее простом, с точки зрения машины, - методе неопределенных коэффициентов. Перед тем как приступить к практике, просмотрим "по диагонали" немного теории, дабы освежить в памяти предмет обсуждения.

Итак, рациональной дробью назовем отношение двух алгебраических многочленов с вещественными коэффициентами:


Дробь называется правильной, если степень P(x) меньше степени Q(x), и неправильной в противном случае. Простейшей называется правильная дробь, знаменатель которой представляет собой неприводимый (значит не имеющий корней) над некоторым полем (в нашем случае — поле действительных чисел) многочлен.

Для простых (правильных) дробей с действительными коэффициентами справедлива следующая теорема о разложении на сумму простейших:

Пусть (1) — правильная рациональная дробь с действительными коэффициентами, знаменатель которой имеет вид:



тогда для этой дроби справедливо следующее разложение на сумму простейших дробей:


где индексированные переменные B, M, N — некоторые вещественные постоянные (может быть, равные нулю).

Для определения конкретных значений сих коэффициентов следует привести равенство к общему знаменателю и сравнить коэффициенты при одинаковых степенях x в числителе. Т.е. по сути дела решить систему линейных уравнений. Используется эта конструкция по большей части при вычислении интегралов, т.к. таким образом интеграл произвольной рациональной функции сводится, по сути дела, к сумме табличных интегралов. Этакая лекция по мат. анализу получилась.

Сразу скажу тем, кому вообще лень что-то делать по этому поводу. Maple делает все, что мы сейчас напишем, одной операцией:

> сonvert(rfun, parfrac, x);

И все. Спросите: зачем этот велосипед? Цель — не конечный результат, а идея и методы ее реализации на Maple. Гораздо интереснее получается посмотреть на целую программу, реально работающий универсальный алгоритм, делающий конкретно нечто, чем просто читать обрывки help-ов под каждую команду языка на английском, не понимая в принципе, как это все связать воедино. Ясное дело, профессионалу, прочувствовавшему Maple, будет неинтересно читать подробные объяснения по поводу использованных функций языка, однако для изучающего систему “не совсем новичка”-математика это будет крайне полезно. Постараюсь в процессе показать читателю свое разумение философии пакета.

Как всегда первый вопрос: с чем работаем? Действительно, для отладки алгоритма необходимо создать хоть несколько рациональных дробей. Руками писать неудобно, поэтому даже этот этап “сгрузим” на машину. Итак, генерируем произвольную рациональную дробь (необязательно правильную) вида:



где s i , t i случайные числа от 1 до 3, P 21 (x), q 1i (x), q 2i (x) — случайные полиномы степени 21, 1, 2 соответственно, причем q 2i (x) не имеют решений над полем действительных чисел.

> restart:
> readlib(randomize):
(а) > randomize():
(б) > d1:= rand(1..3):
> d2:= rand(2..7):
(в) > px:= randpoly(x, degree=21, coeffs=rand(-7..7), terms = 9):
(г) > for i from 1 to 3 do
> q[i]:= randpoly(x, degree=1, coeffs=rand(-7..7))^d1():
> q[i+3]:= (x^2 + x + d2())^d1():
> od:
(д) > rfun:= px/product(q[k], k=1..6);

Разберемся, что тут мы с вами наворотили. Итак, сначала подробно остановимся на генерации случайных целых чисел в системе Maple. (а) — здесь мы заставляем генератор случайных чисел привязаться к системному времени. Если этого не сделать, то генерируемая последовательность будет каждый раз одинаковой. Вызов просто функции rand() без аргументов возвратит двенадцатизначное случайное натуральное число. В большинстве случаев это ну совсем неудобно. Можно это дело исправить, передавая функции один аргумент: rand(n) , что приведет к генерации числа из полуинтервала [0, n) . Зачастую и этого недостаточно для решения поставленной задачи. Можно еще более сузить “область значений” — (б) . Только в этом случае в d1 вернется отнюдь не число, а ссылка на процедуру, вызов которой приведет к генерации случайного числа из заданного отрезка. Произвольный полином максимальной степени 21 степени с коэффициентами из отрезка [-7,7] и девятью членами получим в (в) . Дальше интереснее — нужно изготовить знаменатель. “Сделаем” его в виде произведения трех многочленов первой, и трех — второй степени. Причем по определению многочлены второй степени не должны иметь действительных корней. Реализующая эту задачу конструкция (г) очевидна и в пояснениях не нуждается. И наконец, собрав числитель и знаменатель в одно целое, в (д) получим нашу рациональную дробь. Выражение product(q[k], k=1..6); является формальным переводом на язык Maple записи:


Знаменатель полученной функции таким образом уже будет разложен на множители. Для того чтобы задача приняла более правдоподобный характер, вместо (д) можно реализовать следующее:

> rfun:= px/expand(product(q[k], k=1..6), x);

Expand раскроет скобки и приведет подобные слагаемые относительно переменной x в знаменателе. В итоге получится настоящая рациональная дробь. То, что коэффициенты целые, — общности задачи не ограничивает — с таким же успехом мы могли сгенерировать их и иррациональными.

Если запустить все написанное, исключив строку (а) (генератор случайных чисел будет стандартно инициализирован), получится в точности:


Заметили? В знаменателе появилась “не заказанная” шестая степень. И вместо шести множителей получилось только пять. Ну и что? Просто два “произвольных” многочлена полностью совпали (и по степеням тоже). На что только ни способен генератор случайных чисел в Maple! Результат раскрытия можно посмотреть на рисунке — там он выглядит куда меньше.

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


где Z(x) — целая часть, а R(x) не делится на Q(x) . Сделаем это следующим образом:

> fracpart:= rem(numer(rfun), denom(rfun), x, 'zpart');

Заведем переменную fracpart и zpart соответственно для дробной и для целой части рациональной дроби. Процедура-функция rem возвращает остаток от деления многочленов как основной результат. Третий (необязательный) параметр — имя переменной, “в которую будет вычислена” целая часть. Совершенно аналогично действует функция quo , где основным результатом является целая часть от деления. Здесь функции numer и denom соответственно дают доступ к числителю и знаменателю дроби.

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

> denomx:= factor(denom(rfun));

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



Разделим текущую подзадачу на два этапа: “изготовление” списка знаменателей будущих простейших дробей и запись самого разложения. Для реализации первого этапа нам понадобится написать процедуру-функцию, которая бы занималась преобразованием выражения вида A n в упорядоченный список вида [A, A 2 , A 3 , ..., A n ] .

(а) > transpol:= proc (p: polynom) local j;
(б) > if degree(p, x) <= 1 then
> p;
> else
(в) > if not type(op(2, p), numeric) then
> p;
> else
(г) > seq(op(1, p)^p, j=1.. op(2, p));
> fi;
> fi;
> end:

В (а) объявим имя функции, тип и количество передаваемых параметров, а также локальные переменные в поле local . Результатом работы функции будет результат последней выполненной операции. Теперь опишем сам алгоритм. Если была передана константа либо многочлен первой степени, то вернется он же — (б) . Дальше получим и проанализируем тип объекта op(2, p) . Здесь я обращаюсь к многочлену p как к списку. Maple позволяет работать почти с любым из своих объектов как со списком. После проверки (б) у нас останется лишь три варианта для p : (x 2 +bx+c), (x 2 +bx+c) n , (ax+b) n . Их op(2, p) будет соответственно равен x 2 , n, n . В первом случае (наш “op” — не число) придется возвратить параметр в первозданном виде — это просто квадратный неприводимый трехчлен, а в остальных — осуществляем разложение (г) — формируем нужную последовательность. Далее приготовим список знаменателей будущих простейших дробей, используя только что написанное:

> ds:= [seq(transpol(op(k, denomx)), k=1.. nops(denomx))];

В нашем конкретном случае результат с точностью до расположения элементов списка будет выглядеть следующим образом:

ds:= [2, 2x+3, 5x-4, 4x-1, (4x-1) 2 , (4x-1) 3 , x 2 +x+2, x 2 +x+4, (x 2 +x+4) 2 , ...,(x 2 +x+4) 6 ]

Записать разложение с неопределенными коэффициентами, имея такую прелесть, ничего не стоит:

> rxn:= 0:
> lastvar:= 1:
> for i from 1 to nops(ds) do
(а) > if degree(op(1, op(i, ds)), x) = 1 then
(б) > rxn:= rxn + (A[lastvar])/op(i, ds);
> else
(в) > rxn:= rxn + (A[lastvar]*x+A[lastvar + 1])/op(i, ds);
> lastvar:= lastvar + 1;
> fi;
> lastvar:= lastvar + 1;
> od:

Заставив maple “переварить” сию конструкцию в нашем конкретном случае, получим (опять же с точностью до порядка) нечто тождественное следующему:

Я намеренно не стал “наводить красоту” на результат и “выписывать” именно A i , M j , N j как это сделано в определении разложения и записал только лишь A k -ые, так как смысл от этого не меняется, а работы хоть на четыре строчки, но меньше.

Теперь все сначала и по порядку. Заведем переменную rxn , в которую после запишем разложение. Счетчик lastvar уже использованных индексов коэффициентов установим в значение 1 (следующий, не использованный индекс). Далее, пробегая по списку ds знаменателей будущих простейших дробей, анализируем их степень. Собственно сама реализация такого анализа (а) может показаться довольно странной. Со встроенной функцией degree все понятно — она возвратит степень многочлена относительно переменной, переданной в качестве второго параметра. Что же значит запись op(1, op(i, ds))? Так как вариантов здесь только два, то их и рассмотрим. Если op(i, ds) — выражение вида (x 2 +bx+c) n либо (x+d) n , то op(1, op(i, ds)) вернет то, что находится в скобках. В другом случае — x 2 +bx+c либо x+d (скобок нет) — такая композиция возвратит высший член многочлена (он записан в лексикографическом виде). Таким образом реализуется определение степени знаменателя без учета кратности . А дальше, в зависимости от этого формируется числитель степени на единицу меньшей. За что люблю Maple, так это за (б) и (в) . Ну где вы видели, чтоб вот так “на ходу” можно было “собрать” переменную? А здесь возможно и такое. Естественно, использовав очередной индекс, необходимо увеличить значение счетчика.

Итак, нечто весьма похожее на разложение, приведенное в теореме, мы получили. Теперь дело за малым — нужно вычислить эти самые A k -ые. Сделаем это так: приведем полученное разложение к общему знаменателю, разберемся с подобными и соберем коэффициенты перед x i , где i = 0 ... 21 (в нашем случае) в числителе :

> f:= collect(numer(rxn), x):
> for i from 0 to degree(f, x) do
> cundef[i]:= coeff(f, x, i):
> od:

Функция numer , вернув числитель, “по дороге” приведет rxn к общему знаменателю, collect как раз и повыносит за скобки x i . В переменные (не массив!) cundef i выделим с помощью функции coeff (третий параметр — степень переменной, остальные два очевидны) эти самые коэффициенты. Их количество будет равно степени f плюс один (нулевая). Зачем это надо? А что у нас во fracpart? Именно — то же самое, но коэффициенты определенные. Что делаем? Составляем систему линейных уравнений и решаем относительно наших A k -ых. Единственность решения такой системы доказана до нас, посему спокойно пишем дальше:

> b:= collect(fracpart, x):
(а) > for i from 0 to degree(f, x) do
> cdef[i]:= coeff(b, x, i):
> od:

Снова собрали в cdef i коэффициенты при x i , но уже из fracpart (определенные). Внимание на (а) — их должно быть столько же, сколько и в первом наборе, иначе система не получится. Сформируем набор переменных ( A k- ые), относительно которых будем решать нашу систему и ее саму:

> vars:= {seq(A[k], k=1..lastvar-1)}:
> eqns:= {seq(cundef[i]=cdef[i], i=0.. degree(f, x))}:
> assign(solve(eqns, vars));

Последняя строка заставит Maple пошевелить мозгами, решить нашу систему относительно наших переменных. Функция solve требует два параметра: первый — это набор уравнений, второй — набор переменных. Результат работы будет представлен в виде опять же набора равенств. На этом этапе присвоения переменным, относительно которых решалась система, вычисленных значений не происходит. Чтобы это все-таки сделать, воспользуемся функцией assign в качестве параметра, которой передается набор равенств. Таким образом вычислены наши неопределенные A k -ые. Так как rxn через них выражается, то на результат можно посмотреть так (см. рисунок):

> zpart + rxn;

Это и есть разложение нашей rfun на сумму простейших дробей.

> simplify(zpart + rxn — rfun);

осуществит проверку тождественности (функция simplify как можно дальше упростит выражение), возвратив 0. Другого и не должно быть, в противном случае алгоритм сработал некорректно, чего я от него никак не жду.

Приведенный пример лишний раз доказывает, что система Maple никак не является просто символьным калькулятором. В отличие от MathCad, где возможности программирования представлены не на высшем уровне, она позволяет реализовывать самые буйные фантазии по части построения сложнейших алгоритмов.

Акександр Муравский

© компьютерная газета