Zakład Podstaw Elektrotechniki i Informatyki

METODY NUMERYCZNE – PROJEKT

Częś

ć

Temat

1/2

Interpolacja wielomianowa Lagrange’a

Opracowali

Rok / gr. lab.

Data wyk. ćw.

1. Paweł Lewicki

2. Piotr Lembryk

1EF-DI / L08

14.05.2012

Streszczenie

Praca przedstawia omówienie interpolacji wielomianowej, nazywanej interpolacją Lagrange’a

wraz z przedstawieniem jej implementacji w języku Turbo Pascal, której działanie zostało następnie
porównane z funkcją „polyfit” programu GNU Octave . Dobór przykładów nie jest przypadkowy i ma
na celu uwypuklenie cech charakterystycznych obu metod, zamieszczonych kolejno we wnioskach.

1.Interpolacja Lagrange’a

Interpolacja wielomianowa, nazywana też interpolacją Lagrange'a, od nazwiska pioniera badań

nad interpolacją Josepha Lagrange'a, lub po prostu interpolacją jest metodą numeryczną przybliżania
funkcji tzw. wielomianem Lagrange'a stopnia n, przyjmującym w n+1 punktach, zwanych węzłami
interpolacji wartości takie same jak przybliżana funkcja.

Interpolacja jest często stosowana w naukach doświadczalnych, gdzie dysponuje się zazwyczaj

skończoną liczbą danych do określenia zależności między wielkościami.

Zgodnie z twierdzeniem Weierstrassa dowolną funkcję y=f(x) ciągłą na przedziale domkniętym

można dowolnie przybliżyć za pomocą wielomianu odpowiednio wysokiego stopnia.

Interpolacja liniowa

Jest przypadkiem interpolacji wielomianowej dla dwóch punktów pomiarowych

i

, dla

których można utworzyć funkcję liniową, której wykres przechodzi przez punkty

i

.

Ogólna metoda

Metoda interpolacji polega na:

wybraniu

punktów

należących do dziedziny , dla których

znane są wartości

znalezieniu wielomianu

stopnia co najwyżej takiego, że

.

Interpretacja geometryczna – dla danych

punktów na wykresie szuka się wielomianu

stopnia co najwyżej , którego wykres przechodzi przez dane punkty

Znajdowanie odpowiedniego wielomianu

Wielomian przyjmujący zadane wartości w konkretnych punktach można zbudować w ten
sposób:

Dla pierwszego węzła o wartości

znajduje się wielomian, który w tym punkcie

przyjmuje wartość

, a w pozostałych węzłach

wartość zero.

Dla kolejnego węzła znajduje sie podobny wielomian, który w drugim węźle
przyjmuje wartość

, a w pozostałych węzłach

wartość zero.

Dodaje się wartość ostatnio obliczonego wielomianu do wartości poprzedniego

Dla każdego z pozostałych węzłów znajduje się podobny wielomian, za każdym
razem dodając go do wielomianu wynikowego

Wielomian będący sumą wielomianów obliczonych dla poszczególnych węzłów jest
wielomianem interpolującym

Dowód ostatniego punktu i dokładny sposób tworzenia poszczególnych wielomianów opisany
jest poniżej w dowodzie istnienia wielomianu interpolującego będącego podstawą algorytmu
odnajdowania tego wielomianu.

Dowód istnienia wielomianu interpolującego

Niech

będą węzłami interpolacji funkcji takimi, że znane są wartości

Można zdefiniować funkcję:

,

taką, że dla

jest wielomianem stopnia (mianownik jest

liczbą, a licznik iloczynem wyrazów postaci

)

Gdy

i

:

Gdy

i

:

(licznik = 0 ponieważ występuje element

)

Niech

będzie wielomianem stopnia co najwyżej , określonym jako:

Dla

Wszystkie składniki sumy o indeksach różnych od są równe zeru (ponieważ dla

, składnik o indeksie jest równy:

.

A więc

z czego wynika, że

jest wielomianem interpolującym funkcję

w punktach

.

2.Implementacje w języku Turbo Pascal

2.1 Interaktywny program

Poniżej znajduje się kod programu wykorzystującego funkcję „lagrange.pas” daną w
ramach materiałów pomocniczych projektu użytą w prostym programie interaktywnym:

program zastosowanie_funkcji_lagrange_interaktywny;
function Lagrange (n : Integer;
var x,f : array of Extended;
xx : Extended;
var st : Integer) : Extended;
var i,k : Integer;
fx,p : Extended;
begin
if n<0
then st:=1
else begin
st:=0;
if n>0
then begin
i:=-1;
repeat
i:=i+1;
for k:=i+1 to n do
if x[i]=x[k]
then st:=2
until (i=n-1) or (st=2)
end;
if st=0
then begin
fx:=0;
for i:=0 to n do
begin
p:=1;
for k:=0 to n do
if k<>i then p:=p*(xx-x[k])/(x[i]-x[k]);
fx:=fx+f[i]*p
end;
Lagrange:=fx
end
end
end;
var
global_n, global_st: integer;
global_x, global_f: array[1..20] of Extended;

global_xx: Extended;

i: Integer;

begin
write('podaj liczbe wezlow n : '); readln(global_n);
writeln('podaj ',global_n,' par "wezel wartosc_funkcji" oddzielonych enterem');
for i:=1 to global_n do readln(global_x[i], global_f[i]);
writeln;
writeln('teraz dla kazdego danego x podana zostanie wartosc wielomianu interpolacyjnego, x=100 -> KONIEC');
repeat
read(global_xx);
writeln(Lagrange(global_n, global_x, global_f, global_xx, global_st));
until global_xx = 100;
end.

Przykładowe działanie:

podaj liczbe wezlow n : 5
podaj 5 par "wezel wartosc_funkcji" oddzielonych enterem
1 2
3 3.3
4 2.55
5 7.2
7 11

teraz dla kazdego danego x podana zostanie wartosc wielomianu interpolacyjnego,x=100 -> KONIEC
1
2.0000000000000000E+0000
3
3.3000000000000000E+0000
4
2.5500000000000000E+0000
4.5
4.0669921875000000E+0000
5
7.2000000000000000E+0000
5.6
1.2327380266666667E+0001
8
-4.4166666666666667E+0001
10
-6.1435000000000000E+0002
55
-2.3604310085714286E+0007
100
-5.3809269796428571E+0008

2.2 Program „automatyczny”, wypisujący wartości interpolującego wielomianu dla

zadanego przedziału z zadaną rozdzielczością:

program zastosowanie_funkcji_lagrange_automatyczny;
function Lagrange (n : Integer;
var x,f : array of Extended;
xx : Extended;
var st : Integer) : Extended;
var i,k : Integer;
fx,p : Extended;
begin
if n<0
then st:=1
else begin
st:=0;
if n>0
then begin
i:=-1;
repeat
i:=i+1;
for k:=i+1 to n do
if x[i]=x[k]
then st:=2
until (i=n-1) or (st=2)
end;
if st=0
then begin
fx:=0;

for i:=0 to n do
begin
p:=1;
for k:=0 to n do
if k<>i then p:=p*(xx-x[k])/(x[i]-x[k]);
fx:=fx+f[i]*p
end;
Lagrange:=fx
end
end
end;
var
global_x, global_f: array[1..20] of Extended;
global_xx, global_n, global_st, i: integer;
j, rozdzielczosc: extended;
begin
write('podaj liczbe wezlow n : '); readln(global_n);
writeln('podaj ',global_n,' par "wezel wartosc_funkcji" oddzielonych enterem');
for i:=1 to global_n do readln(global_x[i], global_f[i]);
writeln;
j:=global_x[1]-3;
rozdzielczosc:=0.1;
repeat
writeln(j:5:2,' ', Lagrange(global_n, global_x, global_f, j, global_st):8:5);
j:=j + rozdzielczosc;
until j>global_x[global_n]+3 ;
end.

Przykładowe działanie(dla zakresu <pierwszy_węzeł-3;ostatni_węzeł+3) i rozdzielczości 0.1 ustalonych w
kodzie:

podaj liczbe wezlow n : 4
podaj 4 par "wezel wartosc_funkcji" oddzielonych enterem
2 4.54
3 8.543
8 -2
10 15

-1.00 7.43341
-0.90 6.17306
-0.80 5.04376
-0.70 4.03923
-0.60 3.15333
-0.50 2.38002
-0.40 1.71340
-0.30 1.14768
-0.20 0.67721
-0.10 0.29646
0.00 -0.00000
0.10 -0.21744
0.20 -0.36102
0.30 -0.43579
0.40 -0.44666
0.50 -0.39842
0.60 -0.29573
0.70 -0.14315
0.80 0.05492
0.90 0.29418
1.00 0.57045
1.10 0.87969
1.20 1.21798
1.30 1.58151
1.40 1.96661
1.50 2.36973
1.60 2.78745
1.70 3.21645
1.80 3.65356
1.90 4.09572
2.00 4.54000
2.10 4.98359
2.20 5.42382
2.30 5.85811
2.40 6.28403
2.50 6.69927
2.60 7.10164
2.70 7.48908

2.80 7.85964
2.90 8.21151
3.00 8.54300
3.10 8.85253
3.20 9.13867
3.30 9.40009
3.40 9.63559
3.50 9.84410
3.60 10.02467
3.70 10.17648
3.80 10.29882
3.90 10.39111
4.00 10.45291
4.10 10.48389
4.20 10.48383
4.30 10.45265
4.40 10.39040
4.50 10.29724
4.60 10.17347
4.70 10.01949
4.80 9.83584
4.90 9.62319
5.00 9.38232
5.10 9.11414
5.20 8.81968
5.30 8.50010
5.40 8.15669
5.50 7.79084
5.60 7.40408
5.70 6.99807
5.80 6.57459
5.90 6.13553
6.00 5.68291
6.10 5.21890
6.20 4.74576
6.30 4.26588
6.40 3.78179
6.50 3.29613
6.60 2.81167
6.70 2.33130
6.80 1.85803
6.90 1.39501
7.00 0.94550
7.10 0.51289
7.20 0.10068
7.30 -0.28748
7.40 -0.64784
7.50 -0.97651
7.60 -1.26949
7.70 -1.52265
7.80 -1.73172
7.90 -1.89234
8.00 -2.00000
8.10 -2.05007
8.20 -2.03780
8.30 -1.95832
8.40 -1.80663
8.50 -1.57761
8.60 -1.26601
8.70 -0.86646
8.80 -0.37346
8.90 0.21859
9.00 0.91545
9.10 1.72298
9.20 2.64716
9.30 3.69412
9.40 4.87009
9.50 6.18143
9.60 7.63462
9.70 9.23628
9.80 10.99314
9.90 12.91205
10.00 15.00000
10.10 17.26409
10.20 19.71155
10.30 22.34974

10.40 25.18612
10.50 28.22830
10.60 31.48401
10.70 34.96109
10.80 38.66752
10.90 42.61139
11.00 46.80091
11.10 51.24445
11.20 55.95045
11.30 60.92752
11.40 66.18438
11.50 71.72985
11.60 77.57290
11.70 83.72263
11.80 90.18824
11.90 96.97907
12.00 104.10457
12.10 111.57434
12.20 119.39807
12.30 127.58560
12.40 136.14688
12.50 145.09199
12.60 154.43114
12.70 164.17465
12.80 174.33297
12.90 184.91667

3. Zestawienie z działaniem funkcji „polyfit”

3.1 Osiem węzłów należących do wykresu funkcji x^7-41,5x^6+728x^5-

6992,5x^4+39689x^3-133036x^2+243712x-188145

n=7
węzły:

x[0]=3,5 x[1]=3,9 x[2]=4,2

x[3]=4,6

x[4]=5,6

x[5]=6,4 x[6]=7,6

x[7]=8,1

f[0]=15 f[1]=19 f[2]=6,8489 f[3]=4,7458 f[4]=21,503 f[5]=9,0135 f[6]=20,895 f[7]=16,35

wynik lagrange: funkcja wyjściowa będąca wielomianem(z założenia)
wynik polyfit: 0.85307x^7-35.623x^6+628.58x^5-6070.9x^4+34633x^3-116630x^2+214540x-
166230

3.2 Cztery węzły należące do funkcji: x^5-25x^4+250x^3-1249x^2+3115x-3098

n=3
węzły:
x[0]=4

x[1]=4,5

x[2]=5 x[3]=5,5

f[0]=2

f[1]=2,21875 f[2]=2 f[3]=2,28125

funkcja: x^5-25x^4+250x^3-1249x^2+3115x-3098
lagrange: funkcja wyjściowa będąca wielomianem(z założenia)
polyfit: 1,25x^3-17,75x^2+83,5x-128

3.3 Trzy węzły należące do wykresu funkcji 3^(x-4)

n=2
węzły:
x[0]=4 x[1]=5 x[2]=6
f[0]=1 f[1]=3 f[2]=9
funkcja: 3^(x-4)
lagrange:

2x^2-16x+33

polyfit:

2x^2-16x+33

komentarz:

- wykresy interpolacji polyfit i lagrange nakładają się na siebie
- jak wydać wielomian nie przybliża przebiegu danej funkcji dla x<4

4.4 Trzy węzły należące do wykresu funkcji log10(x-3)+2

n=2
węzły:
x[0]=2,999 x[1]=4 x[2]=5
f[0]=-1

f[1]=2,301

f[2]=2,47712

funkcja: log10(x-3)+2
polyfit: -1,3473x^2+12,4269x-26,1505
lagrange: prawie ta sama co polyfit ( z różnica na 5 miejscu miejscu rozwinięcia dziesiętnego)
komentarz: wykresy interpolacji polyfit i lagrange nakładają się na siebie

4. Wnioski końcowe

- funkcja „lagrange” interpoluje wielomianem który dla punktów węzłowych przyjmuje dokładnie
wartości funkcji interpolowanej, jest nawet w stanie odwzorować całą funkcję interpolowaną o ile
jest ona wielomianem stopnia n+1 gdzie n to liczba węzłów
- przetwarzanie z poziomu kodu Pascala jest szybsze od interpretowanego w Octave
- wraz ze wzrostem stopnia interpolowanego wielomianu błąd rośnie w wypadku Octave
- procedury interpolacji charakteryzujące się mniejszym błędem są zwykle obliczeniowo bardziej
kosztowne
- podobnie ceną dokładności interpolacji metody lagrange’a jest jej złożoność obliczeniowa
- w wypadku interpolacji lagrange'a o ile format zmiennych nie ogranicza precyzji obliczeń błąd
dla węzłów jest zerowy
- lepiej było użyć dużych liczb by uwidocznić błędy, lub ukazując ja jako względne
- algorytm z którego korzysta funkcja polyfit to "metoda najmniejszych kwadratów"

- metoda najmniejszych kwadratów w przeciwieństwie do lagrange’a nie wskazuje najmniejszego
stopnia wielomianu potrzebnego by przeciąć wszystkie wezly, gdyz jej założeniem jest możliwie
dokładne ustalenie wielomianu o stopniu s<=n
- dokładność procedury polyfit zależy jednocześnie od rozmieszczenia punktów węzłowych
względem siebie oraz zadanego stopnia wynikowego wielomianu,
- przebieg wielomianu wynikowego w znacznej mierze zależy od rozmieszczenia względem siebie
węzłów
- charakterystyka działania funkcji lagrange’a pozwala jaj na zerowy błąd w węzłach i
odwzorowanie tej samej funkcji jeżeli węzły należą do wielomianu(bo bazuje na twierdzeniu
Weierstrassa)
- jeśli węzły byłyby punktami przegięcia funkcji której wykresem jest wielomian, funkcja polyfit
mogłaby odwzorować jego przebieg z zerowym błędem – tak jak robi to funkcja „lagrange” bez
względu na to czy węzły są jednocześnie punktami przegięcia