Zakład Podstaw Elektrotechniki i Informatyki

METODY NUMERYCZNE – PROJEKT

Część	Temat
1/2	Interpolacja wielomianowa Lagrange’a
Opracowali	Rok / gr. lab.
Paweł Lewicki Piotr Lembryk	1EF-DI / L08

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

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

1. 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

1. Zestawienie z działaniem funkcji „polyfit”

   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

1. 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

1. 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

1. 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