Дисциплина:
Раздел:
Выполнил:

Вычислительная математика

Интерполяция. Метод Ньютона.

матфак, 3 курс Катешов С.

                                             www.domath.ru

>    restart;

>    #Итак, первая процедура Rc, входные параметры - A,N, где
#A-массив узловых точек и значений интерполируемой функции в этих
#точках, N-порядок разделительной разности.

>    Rc:=proc(A,N)
               local i,j,S,P,x;


#Необходимо объявить локальные переменные,
#в начале процедуры переменные (пм) i,S яв-ся нулевыми.


i:=0:
S:=0:

#При помощи условного оператора создадим механизм разделительных разностей
#, то есть, если пн i не больше, заданного N, то
#выполняется действия P и S.  


while (i<=N) do

#здесь функция product - произведение слагаемых j=0..N
#функция normal приводит выражение к нормальному виду
#ф-ция unnaply возвращает искомый результат


                P:=unapply(normal(product((x-A[j+1][1]),j=0..N)
                            /(x-A[i+1,1])),x):
                S:=S+A[i+1,2]/P(A[i+1,1]):
#накапливание суммы S
 
                i:=i+1:
 end do:

        return S:

end proc;
#конец процедуры

Rc := proc (A, N) local i, j, S, P, x; i := 0; S := 0;  while i <= N do P := unapply(normal(product(x-A[j+1][1],j = 0 .. N)/(x-A[i+1,1])),x); S := S+A[i+1,2]/P(A[i+1,1]); i := i+1 end do; return S end ...
Rc := proc (A, N) local i, j, S, P, x; i := 0; S := 0;  while i <= N do P := unapply(normal(product(x-A[j+1][1],j = 0 .. N)/(x-A[i+1,1])),x); S := S+A[i+1,2]/P(A[i+1,1]); i := i+1 end do; return S end ...
Rc := proc (A, N) local i, j, S, P, x; i := 0; S := 0;  while i <= N do P := unapply(normal(product(x-A[j+1][1],j = 0 .. N)/(x-A[i+1,1])),x); S := S+A[i+1,2]/P(A[i+1,1]); i := i+1 end do; return S end ...
Rc := proc (A, N) local i, j, S, P, x; i := 0; S := 0;  while i <= N do P := unapply(normal(product(x-A[j+1][1],j = 0 .. N)/(x-A[i+1,1])),x); S := S+A[i+1,2]/P(A[i+1,1]); i := i+1 end do; return S end ...
Rc := proc (A, N) local i, j, S, P, x; i := 0; S := 0;  while i <= N do P := unapply(normal(product(x-A[j+1][1],j = 0 .. N)/(x-A[i+1,1])),x); S := S+A[i+1,2]/P(A[i+1,1]); i := i+1 end do; return S end ...
Rc := proc (A, N) local i, j, S, P, x; i := 0; S := 0;  while i <= N do P := unapply(normal(product(x-A[j+1][1],j = 0 .. N)/(x-A[i+1,1])),x); S := S+A[i+1,2]/P(A[i+1,1]); i := i+1 end do; return S end ...
Rc := proc (A, N) local i, j, S, P, x; i := 0; S := 0;  while i <= N do P := unapply(normal(product(x-A[j+1][1],j = 0 .. N)/(x-A[i+1,1])),x); S := S+A[i+1,2]/P(A[i+1,1]); i := i+1 end do; return S end ...
Rc := proc (A, N) local i, j, S, P, x; i := 0; S := 0;  while i <= N do P := unapply(normal(product(x-A[j+1][1],j = 0 .. N)/(x-A[i+1,1])),x); S := S+A[i+1,2]/P(A[i+1,1]); i := i+1 end do; return S end ...

>    #далее следует основная процедура newton
#A-входной параметр процедуры

>    newton:=proc(A)
               local N,Nmax,L,x,i:
#локальные пн процедуры


#первоначальное значение пн L инициализируется, приравнивается к   
#значению интерпол. функции в первом узле.


L:=A[1][2]:



#обозначим размерность списка A
Nmax:=`linalg/vectdim`(A):


#ниже заложен цикл, где и происходит основные вычисления  
#процедуры, при использовании разработанной ранее процедуры #Rc
 
for N from 1 to Nmax-1
                     do
                     L:=L+Rc(A,N)*product(x-A[i+1][1],i=0..N-1):
                    
 #внутри цикла происходит накапливание суммы L.
end do:


#и в конце функция unapply приводит выражение L к нормальному виду.
unapply(collect(L,x),x):
end proc;

newton := proc (A) local N, Nmax, L, x, i; L := A[1][2]; Nmax := `linalg/vectdim`(A); for N to Nmax-1 do L := L+Rc(A,N)*product(x-A[i+1][1],i = 0 .. N-1) end do; unapply(collect(L,x),x) end proc
newton := proc (A) local N, Nmax, L, x, i; L := A[1][2]; Nmax := `linalg/vectdim`(A); for N to Nmax-1 do L := L+Rc(A,N)*product(x-A[i+1][1],i = 0 .. N-1) end do; unapply(collect(L,x),x) end proc
newton := proc (A) local N, Nmax, L, x, i; L := A[1][2]; Nmax := `linalg/vectdim`(A); for N to Nmax-1 do L := L+Rc(A,N)*product(x-A[i+1][1],i = 0 .. N-1) end do; unapply(collect(L,x),x) end proc
newton := proc (A) local N, Nmax, L, x, i; L := A[1][2]; Nmax := `linalg/vectdim`(A); for N to Nmax-1 do L := L+Rc(A,N)*product(x-A[i+1][1],i = 0 .. N-1) end do; unapply(collect(L,x),x) end proc
newton := proc (A) local N, Nmax, L, x, i; L := A[1][2]; Nmax := `linalg/vectdim`(A); for N to Nmax-1 do L := L+Rc(A,N)*product(x-A[i+1][1],i = 0 .. N-1) end do; unapply(collect(L,x),x) end proc
newton := proc (A) local N, Nmax, L, x, i; L := A[1][2]; Nmax := `linalg/vectdim`(A); for N to Nmax-1 do L := L+Rc(A,N)*product(x-A[i+1][1],i = 0 .. N-1) end do; unapply(collect(L,x),x) end proc
newton := proc (A) local N, Nmax, L, x, i; L := A[1][2]; Nmax := `linalg/vectdim`(A); for N to Nmax-1 do L := L+Rc(A,N)*product(x-A[i+1][1],i = 0 .. N-1) end do; unapply(collect(L,x),x) end proc


#продемонстрируем использование данной процедуры на примере следующих данных.

>    XY:=[[-3,5],[-2,3],[0,6],[1,4],[2,7],[3,15]];

XY := [[-3, 5], [-2, 3], [0, 6], [1, 4], [2, 7], [3, 15]]

>    newton(XY)(t);

6-29/360*t^5+5/36*t^4+85/72*t^3-29/36*t^2-73/30*t

>    #функция plot представит график
plot(

[newton(XY)(t),
XY],t=-3..4,color=[green,blue], style=[line,point],title="Интерполяционный многочлен Ньютона",legend=["многочлен","Факт. значения"]);



#обратите внимание, что на одном графике отображены и сам многочлен Ньютона, и фактические значения [x,y], что вполне дает наглядное представление полученного результата

>    #теперь можно сохранить процедуры
#и использовать их по назначению в любом приложении Maple

>    save Rc,newton "newton.m";

#функция Save сохранит обе процедуры в модуле newton.m.
#теперь всякий раз, когда захотите использовать данную процедуру
# в своих программах укажите в самом начале -read "newton.m";-,
#то есть подключите данный модуль, модуль должен находится в той же папке, что и сама программа