====== n-мерный полином, способ задания и его производная ======
В программировании довольно часто приходится работать с полиномами.
Будь то простое квадратное уравнение или многочлен //n//-й степени его все-равно
нужно как-то програмно задать. Рассмотрим один из способов задания полиномов:
Общий вид:
**Pn(x)=a[n]*xn+a[n-1]*xn-1+...+a[1]*x+a[0]**
Если вы захотите его использовать в этом виде то ничего хорошего у вас не выйдет т.к. прийдется использовать функцию возведения в степень, что в свою очередь потянет за собой огромное кол-во умножений (**2n-1**), а операция умножения, как известно, далеко не самая быстрая.
\\ Но это еще не все... Как говорит теория ошибок, погрешность результата прямопропорциональна кол-ву сомножителей т.е.
чем больше вы делаете умножений, тем выше погрешность результата.
\\ Что же делать? Как уменьшить число умножений? На самом деле все просто.
Вспомним школьную математику и вынесем //x// за скобки:
**Pn(x)=(...(a[n]*x+a[n-1])*x+...+a[1])*x+a[0]**
Просмотрев это выражение видно, что мы избавились от операции возведения
в степень и сократили число умножений до **n **раз, а следовательно, повыслили точность и скорость вычислений.
\\ Как же его посчитать? Тут тоже все просто: нужно представить коэффициенты полинома в виде массива и организовать примитивный цикл:
p:=0;
p:=p*x+a[n];
p:=p*x+a[n-1];
. . .
p:=p*x+a[1];
p:=p*x+a[0];
Посмотрим как это можно реализовать на Паскале:
const
anMax = 3;
An: array[0..anMax] of double = (-5.372,1.2493,0.559,-0.13);
{ Тут описан следующий полином: Pn(x)=-5.372 + 1.2493*x + 0.559*x^2 - 0.13*x^3 }
function Pn(x: double): double;
var
i: byte;
Res: double;
begin
Res:=0;
for i:=anMax downto 0 do
Res:=Res*x+An[i];
Pn:=Res;
end;
Это конечно все хорошо, скажете вы, но что-то все равно это не "греет".
\\ Посмотрим теперь пример, как можно легко и просто посчитать //n//-ю производную
заданного полинома.
Вспомним правила вычисления производных:
**1-я (xn)' = n*x(n-1)
2-я (xn)'' = n*(n-1)*x(n-2)
3-я (xn)''' = n*(n-1)*(n-2)*x(n-3)
. . .
i-я (xn)(i) = n*(n-1)*(n-2)*...*(n-i+1)*x(n-i)**
Или на конкретном примере:
**Pn(x) = a[0] + a[1]*x + a[2]*x2 + a[3]*x3
Pn'(x) = a[1] + a[2]*2*x + a[3]*3*x2
Pn''(x) = a[2]*2+ a[3]*3*2*x
Pn'''(x) = a[3]*3*2
Pn''''(x) = 0**
Как видно из примера, чем выше степень производной, тем меньше остается сомножителей и тем меньше мы используем коффициентов из массива, сдвигаясь каждый
раз на позицию вправо. Кроме того, не сложно заметить, что ростом степени производной падает степень и увеличивается коэффициен перед //x//.
\\ Если со степенью все просто и понятно с первого взгляда: из начальной степени вычитаем степень производной и если результат меньше нуля, то выбрасываем этот сомножитель вместе
с его коэффициентом.
\\ А как же расчитать коэффициент? Если немного присмотреться, то можно заметить что коэффициет равен частному двух факториалов:
** k = a! / (a-i)!** , где:
**k **- коэффициент
**a **- начальная степень
**i **- степень производной
Таким образом, получаем формулу:
**i-я (xn)(i) = n!/(n-i)! * x(n-i)**
== Реализация на Паскале: ==
const
anMax = 3;
An: array[0..anMax] of double = (-5.372,1.2493,0.559,-0.13);
function Func(x: double; d: byte):double; far;
{ где d - степень производной }
var
Res: double;
i: byte;
function frac(a,b: byte): longint;
var
res: longint;
begin
res:= 1;
while a > b do
begin
res:= res*a;
dec(a);
end;
frac:= res;
end;
begin
Res:=0;
for i:= anMax-d downto 0 do
Res:=res*x+An[i+d]*frac(i+d,i);
Func:= Res;
end;