Версия для печати
Нажмите сюда для просмотра этой темы в оригинальном формате
Форум на Исходниках.RU > Все языки: Статьи, заготовки в FAQ > n-мерный полином и его производная


Автор: e-moe 10.03.05, 09:47
Полином n-й степени, способ задания и его 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 раз, а следовательно, повыслили точ-
ность и скорость вычислений.
Как же его посчитать? Тут тоже все просто: нужно представить коэффициен-
ты полинома в виде массива и организовать примитивный цикл:

<{CODE_COLLAPSE_OFF}><{CODE_WRAP_OFF}>
                    p:=0;
                    p:=p*x+a[n];
                    p:=p*x+a[n-1];
                       . . .
                    p:=p*x+a[1];
                    p:=p*x+a[0];


Посмотрим как это можно реализовать на Паскале:

<{CODE_COLLAPSE_OFF}><{CODE_WRAP_OFF}>
    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)

Реализация на Паскале:

<{CODE_COLLAPSE_OFF}><{CODE_WRAP_OFF}>
    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;


                                (С) Лабинский Николай. ДонНТУ 2005.
                                        I love Borland.


Last one modified at: GMT+2 09:18:23 10.03.2005

Добавлено
Че скажите?

З.Ы. кому не лень, потестите плиз расчет производной ;)

Powered by Invision Power Board (https://www.invisionboard.com)
© Invision Power Services (https://www.invisionpower.com)