METODY NUMERYCZNE

W ELEKTROTECHNICE

2

Metody numeryczne – dział matematyki stosowanej, zajmujący się opracowywaniem i
badaniem metod przybliżonego rozwiązywania problemów obliczeniowych w modelach
matematycznych innych dziedzin nauki

Przykładowe zastosowania:



elektrotechnika – obliczanie parametrów obwodów elektrycznych



medycyna – tomografia komputerowa, opracowywanie nowych leków



chemia – konstruowanie nowych cząsteczek



inżynieria – przemysł samochodowy i lotniczy



informatyka – konstruowanie nowych procesorów

W obrębie klasycznych metod numerycznych możemy wyróżnić m.in. takie zagadnienia jak:



analiza błędów zaokrągleń



interpolacja



aproksymacja



rozwiązywanie równań i układów równań nieliniowych



całkowanie i różniczkowanie numeryczne



rozwiązywanie układów równań liniowych



obliczanie wartości własnych i wektorów własnych macierzy



rozwiązywanie zagadnień dla równań różniczkowych zwyczajnych i cząstkowych



rozwiązywanie równań całkowych i układów równań całkowych

(każdy wynik musi podlegać weryfikacji!)

Reprezentacja liczb w maszynie cyfrowej
Liczby całkowite (stałopozycyjne = stałoprzecinkowe):

i

n

0

i

i

2

e

znak

l











, gdzie e

i

= 0 lub 1

Jeżeli rejestr ma d-bitów, wówczas liczba całkowita n może zawierać się w przedziale

-2

d-1

< n < 2

d-1

-1

przykład: Zapisać liczbę 18 w systemie dwójkowym:

0

1

2

3

4

2

0

2

1

2

0

2

0

2

1

l





















, może być zatem przechowywana w słowie o długości D

+ 1 > 5 bitów jako













a

rozwinieci

cyfr

d

znak

0010010

...

00

0

Liczby rzeczywiste (zmiennopozycyjne):

c

2

m

zn

x







,

gdzie ½ < m < 1, m-mantysa, i-cecha

i

1

i

i

2

e

m















3

W maszynie cyfrowej mantysa zapisywana jest w t-bitach m

t

i

t

1

i

i

t

)

2

t

(

t

2

e

2

e

m























w wyniku zaokrąglenia do t-bitów mantysy

t

t

2

2

1

m

m







Zmiennopozycyjna reprezentacja liczby rzeczywistej x oznaczana jest symbolem rd(x) i jest
równa

t

c

m

2

zn

)

x

(

rd







t

2

x

x

)

x

(

rd







!

Gdy elementy macierzy są rzędu



lub n musimy dokonać przekształcenia całego

układu. Uciekamy się od operacji na małych liczbach.

Nie prowadzi się operacji na małych liczbach, co oznacza, że:











1

x

)

x

(

rd

,

gdzie

t

2







liczbę 2

-t

nazywamy dokładnością maszynową

przykład: Liczba 18,5 daje się przedstawić w postaci:





6

5

4

3

2

1

)

2

1

2

0

2

1

(

2

1

2

0

2

1

2

0

2

0

2

1

2

x

0

1

2















































może być więc przechowywana w słowie o długości d + 1 > 11 bitów o t > 6 bitach mantysy:

















mantysy

bitów

t

cechy

bitów

t

d

c

znak

znak

0

...

1001010

00101

...

000

0

0



0,5781525

W sytuacji, gdy elementy macierzy są < 1, np. mV, musimy wtedy dokonać przeskalowania
całego układu (np. poprzez pomnożenie macierzy przez jakąś liczbę).

Błędy obliczeń
Przy obliczeniach wykonywanych na maszynach cyfrowych spotyka się trzy podstawowe
rodzaje błędów:



błędy wejściowe (danych wejściowych) – występujące, gdy dane wprowadzane do
pamięci maszyny cyfrowej odbiegają od dokładnych wartości tych danych. Źródłem

4

tych błędów jest najczęściej skończona długość słów binarnych reprezentujących
liczby i w związku z tym nieuniknione jest wstępne ich zaokrąglenie



błędy odcięcia – powstają podczas obliczeń na skutek zmniejszenia liczby działań, na
przykład przy obliczaniu sum nieskończonych. Chcąc obliczyć wartość wyrażenia e

x

równego szeregowi:

...

x

!

N

1

...

x

2

1

x

1

N

2











zastępuje je sumą częściową o

odpowiednio dobranej wartości N:

N

x

N

x

x

!

1

...

2

1

1

2









Jeżeli liczba N będzie niedostatecznie duża, to uzyskana w ten sposób wartość liczby e

x

będzie obliczona niedokładnie, a popełniony w ten sposób błąd jest właśnie błędem odcięcia.
Błędy tego typu powstają też często przy obliczaniu wielkości będących granicami. Podobnie
zastąpienie wartości pochodnej funkcji jej ilorazem różnicowym powoduje powstanie błędu
odcięcia. W wielu przypadkach daje się uniknąć błędu wejściowego i odcięcia przez
ograniczenie danych.

Lemat Wilkinsona:
Błędy zaokrągleń powstające podczas wykonywania działań zmiennopozycyjnych są
równoważne zastępczemu zaburzeniu liczb, na których wykonujemy działania. W przypadku
działań arytmetycznych otrzymujemy:

)

1

(

x

)

1

(

x

)

x

x

(

fl

2

2

1

1

2

1















2

1

3

3

2

1

2

1

)

1

(

x

x

)

x

x

(

fl





















1

2

5

2

1

4

5

2

1

2

4

1

2

1

)

1

(

x

x

x

)

1

(

x

x

x

fl









































gdzie

ε

i

są co do modułu niewiększe niż dokładność maszyny

ε

. Przedrostek fl oznacza, że są

to wyniki działań wyznaczone przy użyciu maszyny cyfrowej.

Oszacowanie błędów w obliczeniach iteracyjnych:
Wiele algorytmów obliczeniowych polega na wyznaczeniu ciągu liczb zbieżnego do
poszukiwanej wartości.

Przykład:
wyznaczyć wartość funkcji

x

y



, przekształcając do równania iteracyjnego. Można to

uczynić obliczając ciąg y

0

, y

1

, y

2

..., gdzie y

0

jest dowolnie wybraną liczbą, natomiast każdy

element ciągu y

i

jest dany wzorem:





















1

i

1

i

i

y

x

y

2

1

y

, dla i = 1, 2, 3, ...

(niewiadoma musi być po obu stronach równania)

2

x

y



x

y

2



x

y

y





y

x

1

y



, współczynnik przed funkcją (tutaj 1) musi być < 1 dlatego:

y

x

y

y

2





5

2

:

y

y

x

y

2

























1

i

1

i

i

y

y

x

2

1

y

- iteracja i-ta

np.

































1

i

1

i

i

2

y

y

x

7

2

1

y

y

x

7

)

y

y

2

(

x

7

y

takiego równania nie da się rozwiązać

bo 7/2
ale można zrobić tak:



























1

i

1

i

i

y

7

y

x

7

8

1

y

y

x

7

)

y

7

y

8

(

jeśli założyć, że działania wykonywane są dokładnie, to ciąg y

0

, y

1

, y

2

, ... będzie zbieżny do

liczby x, ponieważ

x

p

y

i

i



, to

x

p

1

p

2

1

y

i

i

1

i



















, a ciąg p

o

> 0



















i

i

1

i

p

1

p

2

1

p

, przy i = 1, 2, 3, ... jest zbieżny do 1

Przykład:
rozwiązanie iteracyjne





































5

4

)

4

5

(

2

3

7

7

8

5

2

4

3

7

x

y

y

y

x

x

y

x

y

x

















































5

1

2

5

4

2

8

1

3

7

7

1

1

i

i

i

i

i

i

y

x

y

y

x

x

wszystkie składniki przy niewiadomych są < 1 (7/8, -7/8, 3/5)

przyjmujemy dowolnie x

i

= 1, y

i

= 2

Działania potrzebne do wyznaczania kolejnych wartości ciągu y

0

, y

1

, y

2

, ... nazywamy

iteracją. Informację o tym czy w kolejnej iteracji nastąpiło przybliżenie do rozwiązania
uzyskujemy obliczając wyrażenie:

x

y

x

y

q

i

i

i







1

W przypadku dokładnych działań, gdy

x

p

y

i

i



i

i

i

i

p

2

1

p

q





, więc q

i

< 0 dla i = 1, 2, 3, ...

jeśli p

i

>1/3. Oznacza to, że maleje wartość bezwzględna różnicy

x

y

1



odpowiednio

dużego numeru iteracji i.

6

Dla obliczeń wykonywanych na maszynie cyfrowej, mamy:



 



































yi

1

1

x

1

y

2

1

)

y

(

fl

i
3

i
2

i

1

i

1

i

, gdzie

3

i
k





, k = 1, 2, 3, ...

Dla liczby

)

y

(

fl

y

1

i

'

1

i







uzyskamy

x

y

x

y

q

i

'

'
i

1

i









Jeżeli y

i

jest już dostatecznie dobrym przybliżeniem liczby

x , a więc p

i

jest bliskie

jedności, to pierwszy składnik powyższego wyrażenia ma pomijalnie małą wartość, natomiast
drugi ze składników wynika z błędów zaokrągleń i może mieć tym większą wartość im
bliższe jedności jest p

i

.

Przykład: dla liczby 4 x=4

000609

,

2

05

,

2

4

05

,

2

2

1

y

05

,

2

5

,

2

4

5

,

2

2

1

y

5

,

2

4

4

4

2

1

y

1

y

3

2

1

0

















































 





Algorytm:

zmienne rzeczywiste x, y, eps
y:= x/4;

{pierwsze przybliżenie}

repeat
y:= 0,5*(y+(x/y));
until abs (x-sqr(y))<=eps

{eps jest zakładaną dokładnością obliczeń)

Uwarunkowanie zadania i stabilność algorytmów.

Załóżmy, że mamy skończoną liczbę danych rzeczywistych, x = ( x

1

, x

2

, ... x

n

), na ich

podstawie chcemy obliczyć skończenie wiele wyników rzeczywistych y = ( y

1

, y

2

, ... y

m

).

Będziemy więc chcieli określić wartość y według odwzorowania y = φ (x), gdzie



: D → R

m

jest ciągłym odwzorowaniem i D



R

m

.

Algorytm to jednoznaczny przepis obliczania wartości odwzorowania φ (x) składający się ze
skończonej liczby kroków.

Jeśli zadanie rozwiązujemy numerycznie to zamiast dokładnymi wartościami x

1

, x

2

, ... x

n

posługujemy się reprezentacjami rd (x

1

), rd (x

2

), ..., rd (x

n

). Operacje elementarne nie są więc

realizowane dokładnie tylko przez odwzorowanie zastępcze fl.

Ta niewielka zmiana danych

powoduje, że różne algorytmy rozwiązania tego samego zadania dają na ogół niejednakowe
wyniki. Dużą rolę odgrywa tu przenoszenie błędów zaokrągleń. Często niewielkie zmiany
danych powodują duże względne zmiany rozwiązania zadania. Zadanie takie nazywamy źle
uwarunkowanym. Wielkości charakteryzujące wpływ zaburzeń danych na zaburzenia
rozwiązania nazywamy wskaźnikami uwarunkowania zadania.

Przykład:























20

1

k

0

1

19

19

20

20

k

x

a

x

a

...

x

a

x

a

)

x

(

w

7

Zerami tego wielomianu będą liczby naturalne 1, 2, ..., 20. Gdy zakłócimy np. a

19

= - 210 i

jego wartość wynosi a

19

(ε) = - (210 + 2

-23

), czyli a

19

(ε) = a

19

(1+ε). Otrzymujemy wtedy

wielomian w

ε

(x) = w (x) – 2

-23

, który posiada już pierwiastki zespolone i np. najbliższe

pierwiastkowi 15 wielomianu w(x) są pierwiastki 13,992358137



2,518830070j.

Metodę badania przenoszenia błędów można rozbudować do analizy różniczkowej błędów
algorytmu.
Niech δ

x1

= rd (x

i

) – x

i

dla i = 1, 2, ..., n

 



  

i

n

1

i

i

j

yi

x

x

)

x

(

x

x

rd

























dla j = 1, 2, ..., m

We wzorze tym, czynnikiem określającym wrażliwość y

j

na bezwzględną zmianę δx

i

jest

i

j

x

)

x

(







Analogiczny wzór można wyprowadzić dla przenoszenia się błędów względnych. Jeśli

0

y

j



dla j = 0,1, 2, ..., m i

0

x

i



dla i = 1, 2, ..., n to

i

j

x

i

j

n

1

i

j

j

y

x

)

x

(

)

x

(

x





















, gdzie

i

j

j

j

x

)

x

(

)

x

(

x









nazywamy wskaźnikiem uwarunkowania.

Jeśli wskaźniki uwarunkowania co do wartości bezwzględnej są duże to zadanie jest źle
uwarunkowane.

Przykład:
Badamy uwarunkowanie sumy y = a + b + c





R

c

b

a



,

,



 

 



c

b

a

c

b

a

c

b

a



























3

2

1

1

1

1

=















































c

b

a

c

c

b

a

b

c

b

a

a

c

b

a

c

b

a

3

2

1

, gdzie ε jest

największą z wartości ε

1

, ε

2

i ε

3

, gdzie

c

b

a

c

c

b

a

b

c

b

a

a

















jest wskaźnikiem

uwarunkowania zadania.


Stabilność numeryczna algorytmów

Algorytm jest stabilny, jeżeli posiada tę własność, że małe błędy powstałe w pewnym

kroku algorytmu nie są powiększane w innych krokach oraz nie powodują poważnego
ograniczenia dokładności wszystkich obliczeń. Oznacza to, że algorytm jest numerycznie
stabilny, gdy zwiększając dokładność obliczeń można wyznaczyć (z dowolną dokładnością)
istniejące rozwiązanie zadania. Numeryczną stabilność zadania łatwo sprawdzić rozwiązując
zadanie raz z dokładnością np. 10

-6

, a potem z dokładnością 10

-12

.

8

INTERPOLACJA:
Dana jest funkcja y = f (x),





n

x

x

x

,

0



, dla której znana jest tablica wartości w punktach

zwanych węzłami interpolacji. Należy wyznaczyć taką funkcję W(x), aby:
W(x

0

) = Y

0

, W(x

1

) = Y

1

, ..., W(x

n

) = Y

n

interpolacja funkcji f(x)

Zadaniem interpolacji jest wyznaczenie przybliżonych wartości funkcji zwanej funkcją
interpolową w punktach nie będących węzłami interpolacji. Przybliżoną wartość funkcji
obliczamy za pomocą funkcji zwanej funkcją interpolującą, która w węzłach ma te same
wartości co funkcja interpolowana.
Funkcja interpolująca jest funkcją pewnej klasy. Najczęściej będzie to wielomian
algebraiczny, wielomian trygonometryczny, funkcja wymierna i funkcja sklejana.
Interpolację stosuje się najczęściej, gdy nie znamy analitycznej postaci funkcji (jest ona tylko
stablicowana) lub gdy jej postać analityczna jest zbyt skomplikowana.
Najczęściej stosowaną metodą wyznaczania funkcji W(x) jest jej dobór w postaci kombinacji

liniowej n+1 funkcji bazowych

)

x

(

),...,

x

(

),

x

(

),

x

(

n

2

1

0

















n

i

i

i

x

a

x

W

0

)

(

)

(

wyrażenie to nazywamy wielomianem uogólnionym.
Wprowadzając macierz bazową





)

x

(

),...,

x

(

),

x

(

n

1

0











i macierz współczynników



























n

1

0

a

a

a

A



mamy

A

x

x

W







)

(

)

(

warunek W(x

0

) = Y

0

, W(x

1

) = Y

1

, ..., W(x

n

) = Y

n

można zapisać w postaci układu równań liniowych X · A = Y, gdzie

9













































)

x

(

)

x

(

)

x

(

)

x

(

)

x

(

)

x

(

)

x

(

)

x

(

)

x

(

X

n

n

n

1

n

0

1

n

1

1

1

0

0

n

0

1

0

0









































n

1

0

Y

Y

Y

Y



Jeśli macierz X nie jest osobliwa (da się odwrócić), to:

A = X

-1

·Y

co ostatecznie daje

W(x) = Φ(x) · X

-1

· Y

Interpolacja wielomianowa
W praktyce często używa się bazy złożonej z jednomianów

,

)

(

,

,

)

(

,

1

)

(

1

0

n

n

x

x

x

x

x















Baza dla funkcji ciągłych na odcinku skończonym [x

o

, x

n

] jest bazą zamkniętą, tzn., że każda

funkcja tej klasy może być przedstawiona w postaci szeregu złożonego z funkcji bazowych.
Wielomian interpolacyjny ma w tym przypadku postać:

n

n

x

a

x

a

x

a

a

x

W













2

2

1

0

)

(

dodatkowo musi być spełniony warunek:

n

n

n

n

n

1

0

1

n

1

n

1

1

0

0

n

0

n

0

1

0

y

x

a

x

a

a

y

x

a

x

a

a

y

x

a

x

a

a











































Powyższy układ równań ma jedyne rozwiązanie względem a

1

, jeżeli wartości x

0

, x

1

, ..., x

n

są

od siebie różne

n

n

n

n

n

x

x

x

x

x

x

X















1

1

1

det

1

1

0

0



Macierz X

-1

dla bazy wielomianowej

,

x

)

x

(

,

,

x

)

x

(

,

1

)

x

(

n

n

1

0















nazywana jest

macierzą Lagrange’a
Zauważyć należy, że każdy zbiór węzłów równoodległych x

i+1

– x

i

= h = const. można

sprowadzić do zbioru podstawowego podstawiając

h

x

x

q

0





, wówczas

 

n

q

,

0



, a macierz

Φ przyjmuje postać





n

2

q

,

,

q

,

q

,

1







10

Interpolacja Lagrange’a
Przedstawiony powyżej sposób podejścia do interpolacji nie jest zbyt efektywny, ponieważ
macierz X jest macierzą pełną i nie zawsze dobrze uwarunkowaną, co oznacza, że
numeryczna procedura jej odwracania może być obarczona dużym błędem.

W interpolacji wielomianowej Lagrange’a dla n+1 węzłów interpolacji

)

y

,

x

(

,

),

y

,

x

(

,

),

y

,

x

(

),

y

,

x

(

n

n

i

i

1

1

0

0





przyjmuje się funkcje bazowe w postaci

)

x

x

(

)

x

x

)(

x

x

)(

x

x

)(

x

x

(

)

x

(

)

x

x

(

)

x

x

)(

x

x

(

)

x

x

)(

x

x

)(

x

x

)(

x

x

(

)

x

(

)

x

x

(

)

x

x

)(

x

x

)(

x

x

(

)

x

(

)

x

x

(

)

x

x

)(

x

x

)(

x

x

(

)

x

(

1

n

3

2

1

0

n

n

1

i

1

i

3

2

1

0

i

n

3

2

0

1

n

3

2

1

0













































































Funkcje te są wielomianami stopnia n zbudowanymi w ten sposób, że w funkcji bazowej φ

1

brakuje czynnika (x-x

i

). Zatem wielomian interpolacji wyraża się następującym wzorem:

)

x

x

(

)

x

x

)(

x

x

(

a

...

)

x

x

(

)

x

x

)(

x

x

(

a

)

x

x

(

)

x

x

)(

x

x

(

a

)

x

(

W

1

n

1

0

n

n

2

0

1

n

2

1

0



































współczynniki a

0

... a

n

tego wielomianu wyznaczamy z równania:

X · A = Y, przy czym macierz X ma postać:



 





 





 











































1

n

n

0

n

n

1

0

1

n

0

1

0

x

x

...

x

x

0

0

0

x

x

...

x

x

0

0

0

x

x

...

x

x

X















Macierz posiada tylko główną przekątną niezerową w związku z tym układ równań X · A = Y
rozwiązuje się natychmiastowo

)

x

x

(

)

x

x

)(

x

x

(

y

a

)

x

x

(

)

x

x

)(

x

x

(

y

a

)

x

x

(

)

x

x

)(

x

x

(

y

a

1

n

n

1

n

0

n

n

n

n

1

2

1

0

1

1

1

n

0

2

0

1

0

0

0



































11

Można więc wielomian interpolacyjny Lagrange’a zapisać w postaci ułamka:



































n

1

i

n

i

1

i

i

1

i

i

1

i

0

i

n

1

i

1

i

1

0

i

)

x

x

(

)

x

x

)(

x

x

(

)

x

x

)(

x

x

(

)

x

x

(

)

x

x

)(

x

x

(

)

x

x

)(

x

x

(

y

)

x

(

W









lub krócej

 











n

1

i

i

j

j

i

j

i

)

x

x

(

)

x

x

(

y

)

x

(

W

, j = 0, 1, ..., n

Przykład:
Wyznaczyć wielomian interpolacyjny Lagrange’a funkcji f (x) = e

x

w przedziale [0,2 ; 0,5]

mając dane:
f (0,2) = 1,2214, f (0,4) = 1,4918, f (0,5) = 1,6487

009

,

1

91

,

0

724

,

0

1

,

0

3

,

0

)

4

,

0

)(

2

,

0

(

6487

,

1

)

1

,

0

(

2

,

0

)

5

,

0

)(

2

,

0

(

4918

,

1

)

3

,

0

(

2

,

0

)

5

,

0

)(

4

,

0

(

2214

,

1

)

(

2









































x

x

x

x

x

x

x

x

x

W

algorytm do wyznaczania wielomianu Lagrange’a dla podanego punktu:

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;
end;

Należy pamiętać, ze przed przystąpieniem do obliczeń należy sprawdzić czy w wektorze x
wartości się nie powtarzają, ponieważ grozi to wystąpieniem błędu dzielenia przez zero.

Różnice skończone

Dla funkcji stabelaryzowanej przy stałym kroku h = x

i+1

– x

i

można wyprowadzić pojęcie

różnicy skończonej rzędu k

 





 

 









































































k

0

j

1

k

i

k

j

j

i

1

k

1

i

1

k

i

1

k

i

k

i

1

i

2

i

i

1

i

i

i

2

i

1

i

i

y

1

y

y

y

y

...

y

y

2

y

y

y

y

y

y

y

y

12

Przykład: Dla funkcji danej w postaci tablicy zbudować tablicę różnic skończonych.

NR

x

y

Δy

Δ

2

y

Δ

3

y

Δ

4

y

Δ

5

y

0

1,2

1,728

0,469

0,078

0,006

0

0

1

1,3

2,197

0,547

0,084

0,006

0

2

1,4

2,744

0,631

0,090

0,006

3

1,5

3,375

0,721

0,096

4

1,6

4,096

0,817

5

1,7

4,913

Wzory interpolacyjne dla argumentów równoodległych

Dla zbioru węzłów tworzących ciąg arytmetyczny x

0

, x

1

= x

0

+ h, x

2

= x

0

+ 2h, ..., x

n

= x

0

+ nh

dane są wartości funkcji f(x

0

), f(x

1

), ..., f(x

n

). Szukany wielomian interpolacyjny zapisać

możemy w następujący sposób:

W(x) = a

0

+ a

1

q + a

2

q(q – 1) + a

3

q(q – 1)(q – 2) + ... + a

n

q(q – 1)(q – 2)...(q – n +1),

gdzie

h

x

x

q

0





Dla:
x = x

0

: q = 0

x = x

1

: q = 1

x = x

2

: q = 2

...
x = x

n

: q = n

Funkcje bazowe dla wielomianu W(x) przyjęto w postaci:

Φ

0

(x) = 1

Φ

1

(x) = q

Φ

2

(x) = q(q – 1)

Φ

3

(x) = q(q – 1)(q – 2)

Φ

n

(x) = q(q – 1)(q – 2)...(q – n + 1)

Otrzymujemy następujący układ równań:



































































































































n

3

2

1

0

n

3

2

1

0

y

y

y

y

y

a

a

a

a

a

!

n

)

2

n

)(

1

n

(

n

)

1

n

(

n

n

1

0

6

6

3

1

0

0

2

2

1

0

0

0

1

1

0

0

0

0

1



























z którego wyznacza się wartości nieznanych współczynników a

0

, a

1

, a

2

, ..., a

n

13

!

n

y

a

y

a

!

n

...

a

)

1

n

(

n

na

a

...

...

!

3

y

a

y

a

6

a

6

a

3

a

!

2

y

a

y

a

2

a

2

a

y

a

y

a

a

y

a

0

n

n

n

n

2

1

0

0

3

3

3

3

2

1

0

0

2

2

2

2

1

0

0

1

1

1

0

0

0



























































Ostatecznie otrzymujemy I wzór interpolacyjny Newtona w postaci

0

n

0

2

0

0

y

!

n

)

1

n

q

)...(

1

q

(

q

...

y

!

2

)

1

q

(

q

y

q

y

)

x

(

W

























Przykład: Dla zależności f(T)=T log p znaleźć wielomian interpolacyjny stopnia 3 i obliczyć
odchyłki w węzłach interpolacji

T

0,8

1,0

1,2

1,4

1,6

1,8

T log p

0,3566

0,9383

1,5598

2,2169

2,9059

3,6229

Tablica:

T

T log p

Δ

Δ

2

Δ

3

błąd [%]

0,8

0,3566

0,5817

0,0398

-0,0042

0

1,0

0,9383

0,6215

0,0356

-0,0037

0

1,2

1,5598

0,6571

0,0319

-0,0039

0

1,4

2,2169

0,6890

0,0280

0

1,6

2,9059

0,7170

1,9

1,8

3,6229

4,2

przyjmijmy

4

T

5

2

,

0

8

,

0

T

q

8

,

0

T

0















1,5002

1,791T

0,7225T

0,075T

W(T)

2)

1)(q

0,0007q(q

1)

0,0189q(q

0,5817q

0,3566

W(T)

2

3

























Zauważmy, że I wzór Newtona daje dobrą dokładność w pobliżu punktu T

0

, natomiast dla

punktów leżących niżej w tabeli błąd wzrasta. Jeśli chcielibyśmy wyzerować odchyłki we
wszystkich węzłach tabeli trzeba by podwyższyć rząd interpolacji i w efekcie wzór
empiryczny staje się wtedy praktycznie mało przydatny. Wzory interpolacyjne stosuje się też
do obliczania wartości funkcji w punktach pośrednich tabeli (do zagęszczenia). Aby obliczyć
T log p dla T = 1,3 z dokładnością do Δ

2

wystarczy przyjąć T

0

= 1,2, a wówczas q = 0,5 i

znajdujemy

88437

,

1

0319

,

0

25

,

0

6571

,

0

5

,

0

5598

,

1

)

3

,

1

(

W













.

14

Zadanie takich wartości T = 1,65 byłoby już niewykonalne, ponieważ brakuje różnic. Widać
tu ograniczenie tzw. interpolacji w przód, określonej I wzorem Newtona. Do interpolacji
wstecz służy II wzór Newtona. Szukamy wielomian interpolacyjny

W(x) = a

0

+ a

1

q + a

2

q(q + 1) + a

3

q(q + 1)(q + 2) + ... + a

n

q(q + 1)(q + 2)...(q + n – 1),

gdzie

h

x

x

q

n





Współczynniki a

0

, a

1

, a

2

, ..., a

n

wyznaczamy jak poprzednio i dochodzimy do wzoru:

0

2

2

1

0

!

)

1

)...(

1

(

...

!

2

)

1

(

)

(

y

n

n

q

q

q

y

q

q

y

q

y

x

W

n

n

n





























Przykład: Obliczyć T log p dla T = 1,65 z dokładnością do Δ

2

75

,

0

2

,

0

8

,

1

T

q

8

,

1

T

n















)

65

,

1

(

W

3,6229-0,75*0,717-0,75*0,25*0,028=3,0825

for j:=1 to 2 do
for i:=1 to n-j do
z[i,j]:=z[i+1,j-1]-z[i,j-1];
q0:=(x-x0)/h;
k:=int(q,0);
q:=q0-k;
if k>=n then write (‘Brak danych’)
else wart:=z[k,0]+q*z[k,1]+0.5*q*(q-1)*z[k,2];


Interpolacja trygonometryczna

Załóżmy, że znamy wartości pewnej ciągłej i okresowej funkcji f(x) o okresie 2π w 2n+1
węzłach. Jako bazę interpolacji przyjmujemy zbiór funkcji trygonometrycznych













)

nx

cos(

),

nx

sin(

),...,

x

2

cos(

),

x

2

sin(

),

x

cos(

),

x

sin(

,

2

1

Otrzymujemy zatem wielomian interpolacyjny w postaci

)

nx

cos(

a

)

nx

sin(

b

...

)

x

2

cos(

a

)

x

2

sin(

b

)

x

cos(

a

)

x

sin(

b

2

a

)

x

(

W

n

n

2

2

1

1

0

















zawierający 2n+1 nieznanych parametrów.

Ze względu na uproszczenie obliczeń najistotniejszy jest przypadek interpolacji funkcji
określonej na zbiorze równoodległych węzłów x

i



[0,2π] dobranych według następującej

zależności

15

1

n

2

i

2

x

i







, gdzie i = 0, 1, ..., 2n

Czyli

1

n

2

n

4

x

...,

,

1

n

2

2

x

,

0

x

n

1

0















Warunek interpolacji prowadzi do układu równań liniowych w postaci

 

 

 

 

 

 

 

 

















































































































n

1

0

n

1

0

1

1

2n

2n

1

1

1

1

y

...

y

y

a

...

b

a

nx

cos

nx

sin

...

x

cos

x

sin

2

1

...

...

...

...

...

...

nx

cos

nx

sin

...

x

cos

x

sin

2

1

1

0

...

1

0

2

1

Współczynniki pierwszego wiersza macierzy X wynikają z wartości funkcji sin(hx) i cos (hx)
dla x

0

= 0. Przedstawiony układ równań rozwiązuje się natychmiastowo, ponieważ macierz X

-

1

można wyznaczyć z zależności

T

1

X

1

n

2

2

X







Przykład: Daną funkcję f(x) = 7x – x

2

przybliżyć wielomianem trygonometrycznym

przyjmując n = 3. Współrzędne węzłów interpolacji obliczamy ze wzoru:

1

2n

2iπ

x

i





x

0

0,898

1,795

2,693

3,590

4,488

5,386

y

0

5,478

9,344 11,598 12,242 11,274 8,695

Baza interpolacji – zbiór funkcji:













)

x

3

cos(

),

x

3

sin(

),

x

2

cos(

),

x

2

sin(

),

x

cos(

),

x

sin(

,

2

1

16

Tworzymy macierz X:























































































901

,

0

434

,

0

223

,

0

975

,

0

623

,

0

782

,

0

707

,

0

623

,

0

782

,

0

901

,

0

434

,

0

223

,

0

975

,

0

707

,

0

223

,

0

975

,

0

623

,

0

782

,

0

901

,

0

434

,

0

707

,

0

223

,

0

975

,

0

623

,

0

782

,

0

901

,

0

434

,

0

707

,

0

623

,

0

782

,

0

901

,

0

434

,

0

223

,

0

975

,

0

707

,

0

901

,

0

434

,

0

223

,

0

975

,

0

623

,

0

782

,

0

707

,

0

1

0

1

0

1

0

707

,

0

Elementy macierzy X mnożymy przez 2/7, otrzymaną macierz transponujemy i obliczamy
współczynniki wzoru interpolacyjnego





T

491

,

1

,

147

,

0

,

961

,

1

,

513

,

0

,

923

,

4

,

336

,

1

,

845

,

11

A















17

Aproksymacja

jest to przybliżanie funkcji f(x) zwanej aproksymowaną inną funkcją Q(x) zwaną funkcją
aproksymującą. Aproksymacja bardzo często występuje w dwóch przypadkach:



gdy funkcja aproksymowana jest przedstawiona w postaci tablicy wartości i
poszukujemy dla niej odpowiedniej funkcji ciągłej



gdy funkcję o dosyć skomplikowanym zapisie analitycznym chcemy przedstawić w
„prostszej” postaci

Dokonując aproksymacji funkcji f(x) musimy rozwiązać dwa ważne problemy:



dobór odpowiedniej funkcji aproksymującej Q(x)



określenie dokładności dokonanej aproksymacji

Dobór odpowiedniej funkcji aproksymującej Q(x)
Najczęściej stosowane funkcje aproksymujące są dobierane w postaci wielomianów
uogólnionych będących kombinacją liniową funkcji q(x)









m

0

j

j

j

m

)

x

(

a

)

x

(

Q

Jako funkcje bazowe stosowane są:



jednomiany



funkcje trygonometryczne



wielomiany ortogonalne

Przyjęcie odpowiednich funkcji bazowych powoduje, że aby wyznaczyć funkcję
aproksymującą należy wyznaczyć wartości współczynników, a

0

, a

1

, ..., a

m

.

Określenie dokładności aproksymacji
Aproksymacja funkcji powoduje powstanie błędów i sposób ich oszacowania wpływa na
wybór metody aproksymacji. Jeśli błąd będzie mierzony na dyskretnym zbiorze punktów x

0

,

x

1

, ..., x

n

to jest oto aproksymacja punktowa, jeśli będzie mierzony w przedziale (a,b) –

aproksymacja integralna lub przedziałowa.

Najczęściej stosowanymi miarami błędów aproksymacji są:



dla aproksymacji średniokwadratowej punktowej













n

0

i

2

)

x

(

Q

)

x

(

f

S



dla aproksymacji średniokwadratowej integralnej lub przedziałowej











b

a

2

dx

)

x

(

Q

)

x

(

f

S



dla aproksymacji jednostajnej

)

x

(

Q

)

x

(

f

sup

S

b

,

a

x









We wszystkich tych przypadkach zadanie aproksymacji sprowadza się do takiego
optymalnego doboru funkcji aproksymującej (dla wielomianu uogólnionych zaś do
optymalnego doboru współczynników a

0

, a

1

, ..., a

m

) aby zdefiniowane wyżej błędy były

minimalne.

18

Aproksymacja średniokwadratowa
Niech dana będzie funkcja y=f(x), która w pewnym zbiorze X punktów x

0

, x

1

, ..., x

n

przyjmuje wartości y

0

, y

1

, ..., y

n

. Wartości te znane tylko w przybliżeniu z pewnym błędem

(np. jako wyniki pomiarów). Poszukujemy takiej funkcji Q(x) przybliżającej daną funkcję
f(x), która umożliwi wygładzenie funkcji f(x), czyli pozwoli z zakłóconych błędami danymi
wartości funkcji przybliżonej otrzymać gładką funkcję przybliżającą.
Niech



j

(x), j=0, 1,...,n będzie układem funkcji bazowych. Poszukujemy wielomianu

uogólnionego Q(x) będącego najlepszym przybliżeniem średniokwadratowym funkcji f(x) na
zbiorze X=(x

j

). Funkcja przybliżająca ma postać







m

0

i

i

i

(x)

a

Q(x)



Przy czym współczynniki a

i

są tak określone, aby wyrażenie















n

0

i

2

i

i

i

)

x

(

Q

)

x

(

f

)

x

(

w

)

x

(

Q

)

x

(

f

0

)

x

(

w

i



dla i=0, 1, ..., n

było minimalne. Funkcja w(x) jest z góry ustaloną funkcją wagową.

Aby wyznaczyć współczynniki a

i

oznaczamy odchylenie





































n

0

j

2

j

j

2

m

0

i

i

i

i

j

n

0

j

j

n

1

0

R

)

x

(

w

)

x

(

a

)

x

(

f

)

x

(

w

a

,...,

a

,

a

H



































n

0

j

2

j

j

2

m

0

i

j

i

i

j

n

0

j

j

n

1

0

)R

w(x

)

(x

a

)

f(x

)

w(x

a

,...,

a

,

a

H



gdzie R

j

jest odchyleniem w punkcie x

j

.

Następnie obliczamy pochodne cząstkowe funkcji H względem a

i

. Z warunku

0

a

H

k







,

gdzie k = 0, 1, ..., n

otrzymujemy układ m+1 równań o niewiadomych a

i

zwany układem normalnym

0

)

x

(

)

x

(

a

)

x

(

f

)

x

(

w

2

a

H

j

k

m

0

i

i

i

i

j

n

0

j

j

k





































,

gdzie k = 0, 1, ..., n

0

)

(x

)

(x

a

)

f(x

)

w(x

2

a

H

j

k

m

0

i

j

i

i

j

n

0

j

j

k





































Jeśli wyznacznik tego układu jest różny od zera to rozwiązaniem układu jest minimum funkcji
H. W zapisie macierzowym układ przyjmuje postać

F

D

DA

D

T

T



19

gdzie







































n

m

n

0

1

m

1

0

0

m

0

0

x

...

x

...

...

...

x

...

x

x

...

x

D



























m

1

0

a

...

a

a

A



























)

a

(

f

...

)

x

(

f

)

x

(

f

f

n

1

0

Aproksymacja wielomianowa
Jeżeli za funkcje bazowe przyjmiemy ciąg jednomianów (x

i

), i = 0, 1, ..., n to układ normalny

przyjmuje postać

0

x

)

x

(

a

)

x

(

f

n

0

j

k

j

m

0

i

i

j

i

j

























,

k = 0, 1, ..., n

co po zmianie kolejności sumowania daje



 























n

j

m

i

n

j

k

i

j

i

k

j

j

x

a

x

x

f

0

0

0

)

(

przykład: Odnaleźć zależność między x i y w postaci ax+by=1

x 1 3 4 6 8 9 11 14
y 1 2 4 4 5 7 8

9

Będziemy rozpatrywać odchylenia zarówno wartości x jak i y, ponieważ wiodą one do
różnych wyników. Zakładamy, że dane są obarczone błędami, wówczas minimalizujemy
wyrażenia





2

n

1

i

i

i

y

y











2

n

1

i

i

i

x

x







przy czym

1

y

b

x

a

i

i





























n

1

i

i

1

n

1

i

i

0

y

a

x

na

,

gdzie:

b

a

a

1





,

b

1

a

o











































n

1

i

i

i

1

n

1

i

2

i

0

n

1

i

i

y

x

a

x

a

x

2

n

1

i

i

n

1

i

2

i

n

1

i

n

1

i

i

i

i

n

1

i

n

1

i

2

i

i

0

x

x

n

y

x

x

x

y

a























 

 













2

n

1

i

i

n

1

i

2

i

n

1

i

n

1

i

i

i

n

1

i

i

i

1

x

x

n

x

y

y

x

n

a























 













20

Układ normalny dla danych z tablicy ma rozwiązanie
a

0

= 6/11

i

a

1

= 7/11,

więc ma postać
11y – 7x = 6
Minimalizując w ten sam sposób wyrażenie





2

n

1

i

i

i

x

x







otrzymamy równanie
2x - 3y = -1
Dla obu otrzymanych równań wykresy się pokrywają.

Algorytm

tablice: tab[1..4,1..5], a, p[1..4], s, suma[1..6]
zmienne rzeczywiste: x, y
zmienne indeksowe: i, j, n, k

for j:=1 to n do
begin
write (‘Podaj x i y dla węzła:’);
readln (x:10:4, y:10:4);
s[1]:=x;
for i:=2 to 6 do s[i]:=x*s[i-1];
for i:=1 to 6 do suma[i]:=suma[i]+s[i];
p[1]:=p[1]+y;
for i:=2 to 4 do p[i]:=p[i]+y*s[i-1];
end;
tab[1,1]:=n;
for i:=1 to 4 do
begin
tab[i,5]:=p[i];
for j:=1 to 4 do
begin
k:=i+j;
if not (k=2 then tab[i,j]:=suma[k-2]
end;
end;
RUKL(tab,a,4,5);
writeln (‘Poszukiwane współczynniki to’);
for i:=1 to 4 do writeln (a[i]:12:6);

Aproksymacja za pomocą wielomianów ortogonalnych ze wzrostem stopnia wielomianu,
obliczenia stają się coraz bardziej pracochłonne a ich wyniki niepewne. Problem ten można
usunąć stosując do aproksymacji wielomiany ortogonalne.

Funkcje f(x) i g(x) nazywa się ortogonalnymi na dyskretnym zbiorze punktów x

0

, x

1

, ..., x

n

jeśli







n

0

i

i

i

)

x

(

g

)

x

(

f

21

przy czym











n

0

i

2

0

)

x

(

f











n

0

i

2

0

)

x

(

g

Analogicznie ciąg funkcyjny





),...

x

(

),...,

x

(

),

x

(

)

x

(

m

1

0

m











nazywamy ortogonalnym na zbiorze punktów x

0

, x

1

, ..., x

n

jeśli











n

0

i

i

k

i

j

0

)

x

(

)

x

(

dla j≠k

Zastosowanie tej metody powoduje, że znika jedna z trudności obliczeniowych przy
aproksymacji wielomianowej, mianowicie złe uwarunkowanie macierzy układu normalnego.
Przy aproksymacji wielomianami ortogonalnymi macierz układu normalnego jest macierzą
diagonalną, a jej elementy położone na głównej przekątnej dane są wzorem









n

0

i

i

2

j

jj

)

x

(

a

Załóżmy, że znamy n+1 równoodległych punktów x

i

(x

i

= x

0

+ih, i = 0, 1, ..., n). Za pomocą

przekształcenia liniowego

h

x

x

q

0





przeprowadzimy te punkty w kolejne liczby całkowite od 0 do n poszukujemy ciągu
wielomianów





)

q

(

F

),...,

q

(

F

,

F

)

q

(

F

)

n

(

m

)

n

(

1

)

n

(

0

)

n

(

i



(dolny indeks oznacza stopień i m ≤ n) spełniających warunek ortogonalności

0

)

i

(

F

)

i

(

F

n

0

i

)

n

(

k

)

n

(

j







dla j≠k

przy czym

 

 

 

k

k

2

2

1

1

0

)

n

(

k

q

a

...

q

a

q

a

a

)

q

(

F











gdzie

 

)

1

k

q

)...(

1

q

(

q

q

k









22

Często używa się unormowanego ciągu wielomianów spełniających warunek

1

)

0

(

F

)

n

(

k





,

gdzie k = 0, 1, ..., m

]

k

[

k

]

1

[

1

)

n

(

k

q

b

...

q

b

1

)

q

(

F











Co po przekształceniu daje nam wzór na wielomiany Grama, zwane też wielomianami
Czebyszewa stopni k = 0, 1, ..., m w postaci

  

]

n

[

]

s

[

s

k
s

k
s

k

0

s

s

)

n

(

k

n

q

1

)

q

(

F













,

gdzie k = 0, 1, ..., m

Wzór aproksymacyjny oparty na wielomianach Grama ma postać:























 





m

0

k

m

0

k

0

)

n

(

k

k

)

n

(

k

k

m

h

x

x

F

s

c

)

q

(

F

s

c

y

gdzie,

n

m



oraz











n

0

q

2

)

n

(

k

k

)

q

(

F

s

)

x

(

F

y

c

i

)

n

(

n

0

i

i

k









Przykład:

Na podstawie tablicy wartości f(x) znaleźć najlepszy w sensie aproksymacji
średniokwadratowej wielomian przybliżający tę funkcję.

x

1

1,5

2

2,5

3

y

3

4,75

7

9,75

13

Ponieważ węzły są równoodległe posłużymy się wielomianem Grama. Konstruujemy
wielomian aproksymacyjny dla n = 4. Wartości F

k

(q), q=1, 2,…, 4.

Obliczamy ze wzoru:

]

[

]

[

0

)

(

)

1

(

)

(

n

S

k

S

s

n

k

h

q

s

s

k

s

k

q

F











 























1

)

(

0



q

F

n

q

q

F

2

1

)

(

1





1



n

23

)

1

n

(

n

)

1

q

(

q

6

n

q

6

1

)

q

(

F

2











2



n

)

2

n

)(

1

n

(

n

)

2

q

)(

1

q

(

q

20

)

1

n

(

n

)

1

q

(

q

30

n

q

12

1

)

q

(

F

2





















3



n

)

3

n

)(

2

n

)(

1

n

(

n

)

3

q

)(

2

q

)(

1

q

(

q

70

)

2

n

)(

1

n

(

n

)

2

q

)(

1

q

(

q

140

)

1

n

(

n

)

1

q

(

q

90

n

q

20

1

)

q

(

F

2



































4

n



Zestawienie wyników przedstawia tabela:

q

x

i

y

i

F

0

(q)

F

1

(q)

F

2

(q)

F

3

(q)

F

4

(q)

0

1

3

1

1

1

1

1

1

1,5

4,75

1

0,5

-0,5

-2

-4

2

2

7

1

0

-1

0

6

3

2,5

9,75

1

-0,5

-0,5

2

-4

4

3

13

1

-1

1

-1

1

Następnie korzystając ze wzorów:











n

0

q

2

)

n

(

k

k

)

q

(

F

s

)

x

(

F

y

c

i

)

n

(

n

0

i

i

k









Obliczamy c

i

i s

i

oraz b

i

= c

i

/ s

i

i ze wzoru:







m

j

j

j

m

x

a

x

Q

0

)

(

)

(



otrzymujemy:

c

i

37,5

-12,5

1,75

0

0

s

i

5

2,5

3,5

10

70

b

i

7,5

-5

0,5

0

2

dla

5

,

0

1





x

q

24

Aproksymacja trygonometryczna

Często spotykamy się z przypadkiem, gdy funkcja f(x) jest okresowa. Taką funkcję
wygodniej jest aproksymować, wielomianem trygonometrycznym o bazie:

1, sin x, cos x, sin 2x, cos 2x, …, sin kx, cos kx

Jeżeli f(x) jest funkcją dyskretną określoną w dyskretnym zbiorze równoodległych punktów i
ich liczba jest parzysta i wynosi 2L (dla nieparzystej liczby punktów rozumowanie jest
analogiczne), niech:

L

i

x

i





dla i = 0, 1, …, 2L-1

Baza jest ortogonalna nie tylko na przedziale <0, 2π>, ale też na zbiorze punktów x

i

, przy

czym warunki ortogonalności mają postać:































0

k

m

dla

0

0

k

m

dla

L

k

m

dla

0

)

kx

sin(

)

mx

sin(

1

L

2

0

i

i

i































0

k

m

dla

L

2

0

k

m

dla

L

k

m

dla

0

)

kx

cos(

)

mx

cos(

1

L

2

0

i

i

i

0

)

kx

sin(

)

mx

cos(

1

L

2

0

i

i

i











dla m, k dowolnych, przy czym m i k zmieniają się od 0 do L

Przybliżeniem funkcji f(x) na zbiorze punktów x

i

jest wielomian trygonometryczny:











n

1

j

j

j

0

n

))

jx

sin(

b

)

jx

cos(

a

(

a

2

1

)

x

(

y

, n<L

Współczynniki a

j

i b

j

wielomianu wyznaczamy tak, aby suma kwadratów różnic













1

L

2

0

i

2

i

n

i

)

x

(

y

)

x

(

f

była minimalna.

25

Korzystając z warunku ortogonalności otrzymujemy rozwiązanie układu normalnego

0

)

x

(

)

x

(

a

)

x

(

f

)

x

(

w

2

a

H

j

k

n

0

j

j

i

m

0

i

i

j

j

k







































w postaci:



















1

L

2

0

i

1

L

2

0

i

i

i

i

j

L

ij

cos

)

x

(

f

L

1

jx

cos

)

x

(

f

L

1

a

dla j=1, 2, …, n



















1

L

2

0

i

1

L

2

0

i

i

i

i

j

L

ij

sin

)

x

(

f

L

1

jx

sin

)

x

(

f

L

1

b

Algorytm:

const max: … {zdefiniowanie stałej}
type zakres: 1…max;
wektor: array[zakres] of real;
procedure APR (N, K: zakres; var y, a, b: wektor);
var i, j: integer;
PI, A0, S1, S2: real;
begin
write(‘Podaj wartość y’);
for i:=1 to N do readln(Y[I]);
A0:=0.0;
for j:=1 to N do A0:=A0+Y[I];
A0:=A0/N;
PI:=4.0*arctan(1.0)/N;
for i:=1 to K do
begin
S1:=0.0;
S2:=0.0;
for j:=1 to N do
begin
S1:=S1+Y[J]*cos(PI*J*I);
S2:=S2+Y[J]*sin(PI*J*I);
end;
A[I]:=S1*2/N;
B[I]:=S2*2/N;
end;
end;

26

Szybka transformata Fouriera

Każdą funkcję w tablicy tablicą wartości w n punktach

β

n

2π

x

β



, gdzie



= 0, 1, ..., n-1

można przybliżać wielomianem trygonometrycznym w postaci

 









1

n

0

k

k

jkx

exp

c

)

x

(

f

,

j - czynnik urojony

lub



























m

0

k

1

m

k

k

0

x

)

1

m

(

cos

a

2

1

)

kx

sin(

b

)

kx

cos(

a

2

a

)

x

(

f

gdzie



= 1

2

2

n

m





dla n parzystych



= 1

2

1

n

m





dla n nieparzystych





























1

0

1

0

1

0

)

exp(

)

(

1

1

-

n

...,

1,

0,

k

dla

)

sin(

)

(

2

1

)

cos(

)

(

2

1

n

k

n

k

n

k

jkx

x

f

n

c

kx

x

f

n

b

kx

x

f

n

a



















Aproksymacja za pomocą funkcji sklejonych stopnia 3 (funkcji giętych)

Każdą funkcję

)

(

)

(

3

n

S

x

s





można przedstawić w postaci kombinacyjnej liniowej













1

n

1

i

3
i

i

)

x

(

c

)

x

(

s

,

a ≤ x ≤ b oraz gdzie

3
i



są funkcjami określonymi wzorem:





















































































pozostaych

dla

0

x

,

x

x

dla

)

x

x

(

x

,

x

x

dla

)

x

x

(

3

)

x

x

(

h

3

)

x

x

(

h

3

h

x

,

x

x

dla

)

x

x

(

3

)

x

x

(

h

3

)

x

x

(

h

3

h

x

,

x

x

dla

)

x

x

(

n

1

2

i

1

i

3

2

i

1

i

i

3

1

i

2

1

i

1

i

2

3

i

1

i

3

1

i

2

1

i

1

i

2

3

1

i

2

i

3

2

i

3

3
i

gdzie

ih

a

x

i





,

n

a

b

h





Aproksymacja na przedziale <a,b>

Współczynniki c

i

dobieramy tak, aby odchylenie kwadratowe dane poniższym wzorem było

jak najmniejsze:

27

dx

)

x

(

C

)

x

(

f

I

2

b

a

1

n

1

i

3
i

i





























w celu wyznaczenia minimum funkcji I=I(c

-1

, c

0

, … , c

n+1

) obliczamy pochodne:

j

c

I





j = -1 , 0, … , n+1

i przyrównujemy je do zera. Otrzymane równania tworzą układ n+3 równań z n+3
niewiadomymi c

i

 

















1

n

1

i

b

a

b

a

3

j

3

j

3
i

i

dx

)

x

(

)

x

(

f

dx

)

x

(

)

x

(

C

j= -1, 0, 1, … ,n+1

Funkcje

)

x

(

3
i



są liniowo niezależne na przedziale <a,b>, więc układ równań ma

rozwiązanie w postaci punktu. Funkcji I osiąga minimum. Układ równań możemy zapisać w
postaci:















1

n

1

i

b

a

3

j

i

ij

dx

)

x

(

)

x

(

f

h

1

c

a

j = -1, 0, 1, … , n+1,

gdzie









b

a

3

j

3
i

ij

dx

)

x

(

)

x

(

h

1

a

Aproksymacja funkcji określonej na dyskretnym zbiorze punktów

Niech (x

i

’) dla i= 0, 1, … ,n

1

; n

1

> n + 3 będzie zbiorem punktów, w których dane są wartości

funkcji f(x). Szukane jest minimum wyrażenia:

2

n

0

k

1

n

1

i

k

3
i

i

1

)

'

x

(

c

)

x

(

f

J































W celu wyznaczenia c

i

obliczamy pochodne cząstkowe funkcji J = J (c

-1

, c

0

,… ,c

n+1

)

względem zmiennych c

i

i przyrównujemy je do zera.

Otrzymane równania tworzą układ:

















1

n

1

i

n

0

k

k

3
i

k

i

ij

1

)

'

x

(

)

'

x

(

f

c

b

, gdzie











1

n

0

k

k

3

j

k

3
i

ij

)

'

x

(

)

'

x

(

b

28

Jeżeli wyznacznik tego układu jest różny od zera to rozwiązując go możemy wyznaczyć
współczynniki c

j

takie, że funkcja J osiągnie minimum.

Metody numeryczne rozwiązywania układów algebraicznych równań
liniowych

.

Rozpatrujemy układ n równań liniowych zawierających n niewiadomych

n

n

nn

2

2

n

1

1

n

2

n

n

2

2

22

1

21

1

n

n

1

2

12

1

11

b

x

a

...

x

a

x

a

...

...

...

...

...

...

...

b

x

a

...

x

a

x

a

b

x

a

...

x

a

x

a

























Układ ten można zapisać także w postaci macierzowej A · X = B, gdzie















































































n

2

1

n

2

1

nn

2

n

1

n

n

2

22

21

n

1

12

11

x

...

x

x

X

b

...

b

b

B

a

...

a

a

...

...

...

...

a

...

a

a

a

...

a

a

A

metody dokładne: metoda Cramera, metoda eliminacji Gaussa, metoda Crouta (LU)
metody niedokładne: iteracja prosta, Gaussa-Siedla, metoda sukcesywnej nadrelaksacji (SOR)

O wyborze metody decyduje postać macierzy współczynników A. Jeśli macierz A jest
macierzą pełną to na ogół stosujemy metody dokładne. Jeśli macierz współczynników A jest
macierzą niepełną (znacząca ilość współczynników jest równa zero) wtedy stosujemy metody
iteracyjne.
Macierz A nazywana jest macierzą główną układu, X – wektorem niewiadomych, B –
wektorem wyrazów wolnych. Jeśli macierz główna nie jest osobliwa (det A ≠ 0), to układ
równań jest oznaczony (posiada dokładnie jedno rozwiązanie).
Gdy macierz A posiada niezerowe elementy tylko na głównej przekątnej, to układ równań
rozwiązuje się natychmiastowo.



















n

n

nn

2

2

22

1

1

11

b

x

a

b

x

a

b

x

a





Niewiadome można wtedy obliczyć

ii

i

i

a

b

x



, a

ii

≠ 0, i = 0, 1, ..., n.

29

Łatwo rozwiązuje się trójkątne układy równań



































n

n

nn

2

n

n

2

2

22

1

n

n

1

2

12

1

11

b

x

a

...

...

...

b

x

a

...

x

a

b

x

a

...

x

a

x

a

Z ostatniego równania możemy wyznaczyć x

n

, z przedostatniego x

n-1

nn

n

n

a

b

x



a ogólnie

ii

n

1

i

s

s

is

i

i

a

x

a

b

x











, i = n-1, n-2, ..., 1 przy założeniu, że a

ii

≠ 0, i = 1, 2,

..., n

PRZYKŁAD:



































































6

1

6

x

x

x

2

0

0

1

2

0

1

1

1

3

2

1

zmienną x

3

mamy daną wprost - x

3

= 3

1

a

x

a

x

a

b

x

2,

a

x

a

b

x

11

3

13

2

12

1

1

22

3

23

2

2















ALGORYTM:

zmienne indeksowe

n, i, j, s

zmienne rzeczywiste

suma

tablice

a[1..n,1..n], b[1..n], c[1..n]

for i:=n downto 1 do
begin
suma:=0;
for s:=i+1 to n do
suma:=suma+a[i,s]*x[s];
x[i]:=(b[i]-suma)/a[i,i];
end;

Wzory Cramera
Jesli oznaczymy symbolem W wyznacznik główny układu równań













































n

n

nn

2

2

n

1

1

n

2

n

n

2

2

22

1

21

1

n

n

1

2

12

1

11

b

x

a

...

x

a

x

a

...

...

...

...

...

b

x

a

...

x

a

x

a

b

x

a

...

x

a

x

a

30

nn

2

n

1

n

n

2

22

21

n

1

12

11

a

...

a

a

...

...

...

...

a

...

a

a

a

...

a

a

A

det

W



n

2

n

1

n

2

22

21

1

12

11

j

b

...

a

a

...

...

...

...

b

...

a

a

b

...

a

a

W



to można wykazać, że

W

W

x

j

j



, W ≠ 0, j = 1, 2, ..., n

Metoda ta należy do metod dokładnych. Ze względu na dużą złożoność obliczeniową
praktycznie stosowana do numerycznego rozwiązywania równań.

Metoda eliminacji Gaussa
Metoda polega na zapisaniu układu równań w postaci macierzy C, której n pierwszych
kolumn zawiera elementy a

ij

macierzy głównej A, natomiast kolumnę n+1 tworzą dowolne

wyrazy b

i

.

Zakładając, że C

11



0, odejmujemy pierwsze równanie pomnożone przez

11

1

C

C

i

od i-tego

równania (i=2,3,...,n) a obliczonymi współczynnikami zastępujemy przednie wartości. W
wyniku tej operacji otrzymujemy układ równań i odpowiadającą mu macierz C

1

.

)

1

(

1

,

)

1

(

3

)

1

(

3

2

)

1

(

2

)

1

(

1

,

3

)

1

(

3

3

)

1

(

33

2

)

1

(

32

1

,

2

2

3

23

2

22

1

,

1

1

3

13

2

12

1

11

...

...

...

...











































n

n

n

nn

n

n

n

n

n

n

n

n

n

n

n

C

x

C

x

C

x

C

C

x

C

x

C

x

C

C

x

C

x

C

x

C

C

x

C

x

C

x

C

x

C

Za pomocą wzorów określających nowe współczynniki

j

i

ij

ij

C

C

C

C

C

1

11

1

)

1

(





i=2,3,...,n j=2,3,...,n+1

Jeśli

0

)

1

(

22



C

, to odejmiemy drugie równanie układu pomnożone przez

)

1

(

22

)

1

(

2

C

C

i

od i-tego

równania układu i=3,4,...,n. W wyniku otrzymujemy układ:











































)

2

(

1

,

)

2

(

)

2

(

3

)

2

(

1

,

3

)

2

(

3

)

2

(

33

)

1

(

1

,

2

)

1

(

2

)

1

(

23

)

1

(

22

1

,

1

1

13

12

11

2

...

0

0

:

:

...

:

:

:

...

0

0

...

0

...

n

n

nn

n

n

n

n

n

n

n

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

A nowe współczynniki wyrażone są wzorami

)

1

(

2

)

1

(

22

)

1

(

2

)

1

(

)

2

(

j

i

ij

ij

C

C

C

C

C





i=3, 4, ..., n j=3, 4, ..., n

Kontynuując takie postępowanie po wykonaniu n kroków dochodzimy do układu trójkątnego

31

)

1

(

1

,

)

1

(

)

2

(

1

,

1

)

2

(

3

3

)

2

(

33

)

1

(

1

,

1

)

1

(

2

3

)

1

(

23

2

)

1

(

22

1

,

1

1

3

13

2

12

1

11

...

...

...







































n

n

n

n

n

nn

n

n

n

n

n

n

n

n

n

C

x

C

C

x

C

x

C

C

x

C

x

C

x

C

C

x

C

x

C

x

C

x

C

















































)

1

(

1

,

)

1

(

)

2

(

1

,

3

)

2

(

3

)

2

(

33

)

1

(

1

,

2

)

1

(

2

)

1

(

23

)

1

(

22

1

,

1

1

13

12

11

2

0

0

0

0

0

0

n

n

n

n

nn

n

n

n

n

n

n

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C





















Algorytm dochodzenia od układu trójkątnego do układu równań liniowych











































)

1

(

)

1

(

)

1

(

)

1

(

)

(

,...,

2

,

1

1

,...,

2

,

1

s

sj

s

ss

s

is

s

ij

s

ij

C

C

C

C

C

n

s

s

i

n

s

Przykład



































4

2

3

5

2

3

1

4

3

2

3

2

1

3

2

1

3

2

1

x

x

x

x

x

x

x

x

x



































4

2

3

1

5

1

2

3

1

4

3

2

C

Obliczamy elementy macierzy C

1

2

9

0

2

9

2

13

7

2

13

14

11

31

34

)

1

(

34

13

11

31

33

)

1

(

33

12

11

31

32

)

1

(

32

14

11

21

24

)

1

(

24

13

11

21

23

)

1

(

23

12

11

21

22

)

1

(

22









































C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

32

druki krok s=2

0

13

63

)

1

(

24

)

1

(

22

)

1

(

32

)

1

(

34

)

2

(

34

)

1

(

23

)

1

(

22

)

1

(

32

)

1

(

33

)

2

(

33















C

C

C

C

C

C

C

C

C

C











































0

13

63

0

0

2

13

7

2

13

0

1

4

3

2

2

C

co daje nam macierz C

1

:











































2

9

0

2

9

0

2

13

7

2

13

0

1

4

3

2

1

C

i dalej ze wzorów dla układów trójkątnych mamy x

3

=0, x

2

=-1, x

1

=1

ALGORYTM:

zmienne indeksowe

n, i, j, s

zmienne rzeczywiste

suma

tablice

c[1..n,1..n+1], x[1..n]

for s:=1 to n-1 do
for i:=s+1 to n do
for j:=s+1 to n+1 do
c[i,j]:=c[i,j]-(c[i,s]*c[s,j])/c[s,s];

x[n]:=c[n,n+1]/c[n,n]

for i:=n-1 downto 1 do
begin
suma:=0;
for s:=i+1 to n do
suma:=suma+c[i,s]*x[s];
x[i]:=(c[i,n+1]-suma)/c[i,i];
end;


Metody iteracyjne
Proces iteracyjnego obliczania wrtości x polega na przyjęciu pewnej wartości początkowej
X

o

, zwanej przybliżeniem początkowym, a następnie wykonaniu określonych operacji

arytmetycznych dających w wyniku X

1

(zwanym pierwszym przybliżeniem). Kolejne

przybliżenia oblicza się wykonując te same operacje na ostatnio obliczonym przybliżeniu.

33

Jeżeli ciąg przybliżeń {X

n

} zmierza do X, to proces iteracyjny jest zbieżny i tylko wtedy

metoda iteracyjna jest efektywna. Metody iteracyjnego rozwiązywania układów równań
liniowych można stosować dla układów typu

...













































































nn

n2

n1

2n

22

21

1n

12

11

n

2

1

w

w

w

w

w

w

w

w

w

X

X

X

















Przykład: Sprowadzić układ do postaci AX=B

2x

1

+x

2

+x

3

=4

x

1

=-0,5x

2

-0,5x

3

+2

x

1

+x

2

-x

3

=1



x

2

=-x

1

+x

3

+1

x

1

+x

2

+x

3

=3

x

3

=-x

1

-2x

2

+4

W pierwszej kolejności macierz A zakładamy na dwie macierze D i R przy czym D zawiera
tylko główną przekątną macierzy A, natomiast R pozostałe jej elementy.





































































0

2

1

1

0

1

1

1

0

1

0

0

0

1

0

0

0

2

1

2

1

1

1

1

1

1

2

A

=

D

+

R

Równanie A*X=B sprawdziło się dla

D*X+R*X=B  D*X=-R*X+B

Ostatnie równanie pomnożymy lewostronnie przez D

-1

i otrzymujemy:

E*X=D

-1

*R*X+D

-1

*B

Oznaczając –D

-1

*R=W i D

-1

*B=Z dochodzimy dorównania: X=W*X+Z

Wracając do zadania:











































































































































4

1

4

1

0

0

0

1

0

0

0

2

0

2

1

1

0

1

1

1

0

1

0

0

0

1

0

0

0

2

1

3

2

1

1

3

2

1

x

x

x

x

x

x

ponieważ













































1

0

0

0

1

0

0

0

5

,

0

1

0

0

0

1

0

0

0

2

1

więc jak łatwo sprawdzić mamy taki sam układ jaki znaleziono sposobem „naturalnym”

34

Metoda iteracji prostej (Jakobiego):
Metoda ta dla równania X=W*X+Z polega na przyjęciu zerowego przybliżenia wektora
X=X

o

, a następnie przeprowadzenia obliczeń iteracyjnych za pomocą zależności:

X

i+1

=W*X

i

+Z i=0,1,...

Liczba kroków, które należy wykonać, aby uzyskać wymaganą dokładność rozwiązania, jest z
reguły dość duża i istotnie zależy od przyjęcia punktu startowego X

o

. Punkt ten najczęściej się

dobiera na podstawie dodatkowych informacji o fizycznych aspektach problemu.

Przykład: Metodę iteracji prostej wykorzystano do rozwiązania układu równań będącego
różnicową aproksymacją równania Laplace`a z warunkami brzegowymi I rodzaju. Równanie
to opisuje stacjonarne pole temperatury w pewnym dwuwymiarowym obszarze, na którego
brzegach temperatury są znane. Jako przybliżenie zerowe przyjęto wartość temperatur w
wnętrzu obszaru T=50



C, wymiary prostokątne a=b=0,8 cm.

























































50

0

0

0

0

0

0

0

0

100

50

38

38

38

38

38

25

0

100

63

50

50

50

50

50

38

0

100

63

50

50

50

50

50

38

0

100

63

50

50

50

50

50

38

0

100

63

50

50

50

50

50

38

0

100

63

50

50

50

50

50

38

0

100

75

63

63

63

63

63

50

0

100

100

100

100

100

100

100

100

50

























































0

0

0

0

0

18

14

10

6

0

32

26

19

10

0

42

35

26

14

0

50

42

32

18

0

58

50

39

22

0

68

61

50

31

0

100

94

90

86

82

78

69

50

0

100

100

100

100

100

100

100

100

50

35

Metoda Gaussa-Seidla

Stanowi ulepszenie metody iteracyjnej prostej polegające na wykorzystaniu obliczonych k
pierwszych składowych wektora X

i+1

do obliczenia składowej k+1. Modyfikacja ta znacznie

przyspiesza proces obliczania.

Przykład: Rozwiązując układ równań:









































































































2

1

5

2

0

3

,

0

1

,

0

2

,

0

1

,

0

0

2

,

0

3

,

0

2

,

0

3

,

0

0

1

,

0

2

,

0

1

,

0

2

,

0

0

4

3

2

1

4

3

2

1

x

x

x

x

x

x

x

x

dla uproszczenia zapisu oznaczamy składowe z wektora iteracji „i” symbolem „x”, a
składowe wektora i iteracji i+1 symbolem „u”. W metodzie iteracji prostej układ ten
zapisujemy:

2

3

,

0

1

,

0

2

,

0

1

1

,

0

2

,

0

3

,

0

5

2

,

0

3

,

0

1

,

0

2

2

,

0

1

,

0

2

,

0

3

2

1

4

4

2

1

3

4

3

1

2

4

3

2

1

































x

x

x

u

x

x

x

u

x

x

x

u

x

x

x

u

Przyjmujemy teraz punkt startowy np. X

0

= [3, 6, 4, 3]

T

, składowe tego wektora podstawiamy

do prawych stron równań wyliczając u

1

, ..., u

4

, które w następnej iteracji przejmą rolę punktu

startowego. Kolejne przybliżenia otrzymane tą metodą

































































































































78

,

3

13

,

4

33

,

7

53

,

4

37

,

3

05

,

4

27

,

7

45

,

4

15

,

3

4

08

,

7

4

,

4

2

,

3

4

,

3

1

,

7

2

,

4

3

4

6

3

rozwiązaniem dokładnym jest X = [4,59, 7,49, 4,20, 3,44]

Ten sam układ w metodzie Gaussa – Siedla ma postać:

2

3

,

0

1

,

0

2

,

0

1

1

,

0

2

,

0

3

,

0

5

2

,

0

3

,

0

1

,

0

2

2

,

0

1

,

0

2

,

0

3

2

1

4

4

2

1

3

4

3

1

2

4

3

2

1

































u

u

u

u

x

u

u

u

x

x

u

u

x

x

x

u

36

a kolejne przybliżenia:

































































































































44

,

3

2

,

4

4

,

7

58

,

4

43

,

3

19

,

4

38

,

7

56

,

4

41

,

3

15

,

4

32

,

7

51

,

4

32

,

3

4

22

,

7

2

,

4

3

4

6

3

Formalny zapis algorytmu

Macierz W w równaniu X = W · X + Z przedstawiamy jako sumę dwóch macierzy: górnej
trójkątnej W

u

i dolnej trójkątnej W

l

. Macierz górna trójkątna składa się z elementów macierzy

W leżących nad główną przekątną (pozostałe są zerami), natomiast macierz dolna trójkątna
składa się z elementów macierzy W leżących pod główną przekątną. Macierz W w tym
przykładzie ma zerową główną przekątną

Algorytm metody Gaussa – Siedla realizowany jest następującym wzorem:

Z

X

W

X

W

X

i

i

i

u

i















1

1

czyli mamy:



































































































































2

1

5

2

0

3

,

0

1

,

0

2

,

0

0

0

2

,

0

3

,

0

0

0

0

1

,

0

0

0

0

0

0

0

0

0

1

,

0

0

0

0

2

,

0

3

,

0

0

0

2

,

0

1

,

0

2

,

0

0

4

3

2

1

4

3

2

1

x

x

x

x

u

u

u

u

Metoda sukcesywnej nadrelaksacji (SOR)

Jest ulepszeniem metody Gaussa-Seidla mającym poprawić szybkość zbieżności procesu
iteracyjnego. Istota polega na sukcesywnym wprowadzaniu w miejsce składowych u

i

po

prawej stronie układu wartości

x

i

+ α (u

i

– x

i

)

α – współczynnik nadrelaksacji, α



<1,2>

i – liczba niewiadomych a nie iteracja

Sposób dobory optymalnego współczynnika



nazywanego czynnikiem nadrelaksacyjnym

jest trudny do określenia, chociaż wiadomo, że jest on liczbą z przedziału [1 , 2]. Można
zauważyć, że dla



=1 otrzymuje się metodę Gaussa – Siedla. Według źródeł literaturowych

zaleca się przyjmowanie



rzędu 1,8.

Algorytm metody nadrelaksacyjnej sprowadza się do wzoru:

Z

X

X

X

W

X

W

X

i

i

i

i

i

u

i



















)]

(

[

1

1

1



37

Zbieżność metod iteracyjnych

Jeżeli macierz A kwadratowa oraz x jest wektorem niezerowym, wówczas istnieje pewien
wektor równoległy do Ax, co oznacza, że istnieje pewna stała



taka, że:

Ax =



x

Możemy wtedy zapisać równanie charakterystyczne

(A -



I)x = 0,

gdzie I – macierz jednostkowa

oraz wielomian charakterystyczny ma postać

p(



) = det (A -



I)

Jeżeli wielomian macierzy A jest

n

n



, wektora x jest n, wówczas p jest wielomianem m-

tego stopnia, którego pierwiastki są pojedyncze i zespolone. Jeżeli wektor



jest

pierwiastkiem wielomianu p wówczas:





0

det





I

A



Jeżeli powyższy warunek jest spełniony, oznacza to, że równanie charakterystyczne ma
rozwiązanie x



0.

Zbieżność metod iteracyjnych (2)

Promień spektralny



(A) macierzy A jest zdefiniowany jako:







max

)

A

(

,

gdzie



- oznacza moduł wektora



wówczas równanie charakterystyczne jest zbieżne gdy:

1

)

A

(





Można wykazać, że:

A

macierzy

normę

oznacza

A

gdzie

A

A



)

(



ostatecznie dla metod iteracyjnych można zapisać twierdzenie:


Warunkiem koniecznym i wystarczającym zbieżności procesu iteracyjnego danego wzorem

Z

X

W

X

i

i









1

przy dowolnym wektorze początkowym

)

0

(

X

i dowolnym wektorze Z jest,

aby wszystkie wartości własne macierzy W były co do modułu mniejsze od jedności.

38

Zbieżność metod iteracyjnych

W praktyce jest na ogół wygodniej korzystać z twierdzenia:

Jeżeli norma macierzy W jest mniejsza od jedności to proces iteracyjny

określany wzorem X

i+1

= W · X

i

+ Z jest zbieżny.

Z przytoczonych twierdzeń wynika, że zbieżność procesu iteracyjnego nie zależy ani od
wartości początkowej

)

0

(

X

, ani też od wartości danego wektora Z, a wystarczy jedynie

spełnienie jednego z następujących warunków:







n

j

ij

i

w

1

1

max

norma wierszowa







n

i

ij

j

w

1

1

max

norma kolumnowa









n

i

n

j

ij

w

1

1

2

1

norma energetyczna

Kryterium „stopu” metod iteracyjnych

Obliczenia iteracyjne trwają do momentu kiedy zostanie spełniony warunek

 

 









i

1

i

X

X

gdzie: X

(i+1)

oraz X

(i)

są kolejnymi przybliżeniami rozwiązania X

i+1

= W · X

i

+ Z oraz δ jest

dokładnością obliczeń – dowolna liczba większa od 0















j

1

i

j

i

1

i

x

x

x

,

gdzie ε



<10

-1

, 10

-256

> – przyjęta dokładność obliczeń

39

Numeryczne rozwiązywanie równań nieliniowych

Metoda Newtona

Dany jest układ równań nieliniowych z n niewiadomymi

0

)

x

...

x

(

f

..

..........

..........

0

)

x

...

x

(

f

0

)

x

...

x

(

f

n

1

n

n

1

2

n

1

1







który możemy zapisać w postaci wektorowej

f(x)=0

gdzie:



























n

x

x

x

:

:

1



























)

(

:

:

)

(

)

(

1

x

f

x

f

x

f

n

O funkcjach fi(x

1

...x

n

) dla i = 1, 2, ..., n zakładamy, że mają ciągłe pochodne pierwszego

rzędu w pewnym obszarze zawierającym odosobniony pierwiastek układu równań. Niech

x

(k)

= {x

1

(k)

, ..., x

n

(k)

}

będzie k-tym przybliżeniem pierwiastka a = {a

1

, ..., a

n

} równania.

Dokładną wartość pierwiastka wyraża wzór a = x

(k)

+ ε

(k)

Gdzie

}

,...,

{

)

(

)

(

1

k

n

k









jest błędem pierwiastka przybliżonego

)

(k

x

. Skoro istnieje ciągła

pochodna funkcji f(x) możemy zapisać:

)

(

)

(

)

(

)

(

'

)

(

)

(

k

k

k

x

f

x

f

a

f









Zastępując błąd

)

(k



przyrostem

)

(

)

1

(

)

(

k

k

k

x

x

h







i porównując prawą stronę powyższego

wyrażenia do zera otrzymujemy równanie liniowe:

0

)

)(

(

'

)

(

)

(

)

1

(

)

(

)

(









k

k

k

k

x

x

x

f

x

f

Wzór ten stanowi zapis rekurencyjny dla metody Newtona w postaci macierzowej.

40

Pochodna f’(x) jest macierzą Jakobiego





















































































n

n

2

n

1

n

n

2

2

2

1

2

n

1

2

1

1

1

x

f

...

x

f

x

f

...

...

...

...

x

f

...

x

f

x

f

x

f

...

x

f

x

f

)

x

(

W

)

x

(

'

f

Przyjmując, że

)

(

)

(k

x

W

jest macierzą nieosobliwą możemy równanie liniowe przekształcić do

postaci:

)

(

)

(

)

(

)

(

1

)

(

)

1

(

k

k

k

k

x

f

x

W

x

x











stąd przyjmując dowolną wartość

)

0

(

x

otrzymujemy ciąg kolejnych przybliżeń pierwiastka

równania

)

(

)

2

(

)

1

(

)

0

(

...,

,

,

,

k

x

x

x

x

uzyskanych metodą Newtona.

Metoda iteracji

Dany jest układ równań nieliniowych z n niewiadomymi

)

x

...

x

(

g

x

..

..........

..........

)

x

...

x

(

g

x

)

x

...

x

(

g

x

n

1

n

3

n

1

2

2

n

1

1

1







który możemy zapisać w postaci wektorowej

x=g(x)

gdzie:



























n

x

x

x

:

:

1



























)

(

:

:

)

(

)

(

1

x

g

x

g

x

f

n

Zakładamy, że funkcja g

1

, g

2

, ..., g

n

są rzeczywiste i ciągłe w pewnym otoczeniu

odosobnionego pierwiastka a = {a

1

, ..., a

n

} układu równań.

41

Metoda iteracji polega na stosowaniu następującego wzoru

x

k+1

= g(x

(k)

)

Po przyjęciu przybliżenia

)

0

(

x

, w rezultacie otrzymujemy ciąg kolejnych przybliżeń

...

,

,

,

)

2

(

)

1

(

)

0

(

x

x

x

pierwiastka a równania.

Metoda siecznych

Rozpatrujemy układ równań nieliniowych w postaci



















0

)

x

...,

,

x

,

x

(

f

.

..........

..........

..........

0

)

x

...,

,

x

,

x

(

f

0

)

x

...,

,

x

,

x

(

f

n

2

1

n

n

2

1

2

n

2

1

1

Który możemy ogólnie zapisać w postaci:

f(x)=0

gdzie x jest wektorem n zmiennych.
Metoda iteracji polega na stosowaniu wzoru wyliczającego k-te przybliżenie i-tej zmiennej
tych:

)

(

)

(

)

(

)

(

1

1

1

















k

i

k

i

k

i

k

i

k

i

k

i

k

i

x

f

x

f

x

f

x

x

x

x

Metoda ma dwie wady:
- potrzebne są 2 zestawy punktów startowych: x

0

, f(x

0

), x

1

, f(x

1

)

- wykrywamy jeden wektor x-ów (globalny) (jeden zestaw x-ów), a pierwiastki mogą być
wielokrotne


Przykład: Dany jest układ równań w postaci

















0

2

y

x

0

4

x

2

Napisać program wyznaczający rozwiązanie tego układu.
Niezbędne deklaracje typów:

42

wektor array[1..2] of Extended;

zmienne:

x, X1, Xn, F, F1, Fn: Extended;
iteracja: LongInt;

Funkcja reprezentująca układ równań

function funkcja (z:word):wektor;
begin
funkcja[1]:=sqrt(z[1])-4;
funkcja[2]:=z[1]*z[2]-2;
end;

Realizacja obliczeń:

Procedure oblicz
Begin
x[1]:=4;
x[2]:=7;
X1[1]:=9;
X1[2]:=11;
Repeat
Inc(iteracje);
F:=funkcja(x);
F1:=funkcja(X1);
xn[1]:=x[1]

-

((x[1]-x1[1])*F[1])/F[1]-F1[1];

Xn[2]:=x[2]

-

((x[2]-X1[2])*F[2])/F[2]-F1[2];

Fn:=funkcja(xn);
X1:=x;
x:=Xn;
until (abs(Fn[1]+abs(Fn[2])<=0.0000000001;

Po zakończeniu obliczeń wektor x przechowuje wartości zmiennych, natomiast zmienna
iteracje zawiera liczbę wykonanych iteracji.

43

Rozważania dotyczą numerycznego obliczania wartości:
- całek pojedynczych funkcji ciągłych
- całek pojedynczych funkcji niewłaściwych
- całek funkcji osobliwych
- całek wielokrotnych
- całek Monte Carlo

Całkowanie pojedyncze – metody obliczeń

a) kwadratury z ustalonymi węzłami



kwadratury Newtona-Cotesa



złożone kwadratury Newtona-Cotesa



metody Romberga



metoda adaptacyjna

b) kwadratury Gaussa

- kwadratury Gaussa
- złożone kwadratury Gaussa

Całki pojedyncze właściwe

Całka oznaczona Riemanna
Niech funkcja f(x) będzie określona i ograniczona na przedziale [a,b]. Dokonajmy podziału
przedziału [a,b] takiego, że:

a = x

0

< x

1

< ... < x

n+1

= b

i utwórzmy sumę



  

 

1

c

f

x

x

S

n

0

i

i

i

1

i

n













,

gdzie

1

i

i

i

x

,

x

c





są dowolnymi punktami pośrednimi

Jeżeli ciąg {S

N

} dla





n

jest zbieżny do tej samej granicy S przy każdym przedziale

odcinka [a ; b] takim, że

0

max

1

0











i

i

n

i

x

x

niezależnie od wyboru punktów C

i

to funkcję f(x)

nazywamy całkowalną w przedziale [a ; b], a granicę S ciągu (1) nazywamy:
Całką oznaczoną Riemanna funkcji f o oznaczamy:



b

a

dx

x

f

)

(

Można wykazać, że funkcja ograniczona w przedziale domkniętym jest całkowalna wtedy i
tylko wtedy, gdy jej punkt nieciągłości w przedziale [a ; b] tworzą zbiór miary zero.

Jeżeli przez F(x) oznaczymy funkcję pierwotną funkcji f(x) w przedziale [a,b] (jeżeli F’(x) =
f(x)) to ma miejsce wzór

44

 

 

   

 

2

a

F

b

F

dx

x

f

f

I

b

a









 

 

 

3

f

E

S

f

I

n





, gdzie E(f) – reszta kwadratury

Interpretacją graficzną całki oznaczonej jest pole pomiędzy osią OX, krzywą.
Wzór (1) jest inaczej nazywany kwadraturą. Nazwa pochodzi od sumowania wyrazów ciągu,
które w najprostszym przypadku graficzną interpretację mają w postaci czworokątów.

45

Całka oznaczona Riemanna

Interpretacją graficzną całki oznaczonej jest pole powierzchni pomiędzy osią OX a krzywą.
Rysunek obok to ilustruje

Wyznaczmy S sumę ciągu n-wyrazów, które są określone jako iloczyn wartości funkcji w
punktach x

i



<a,b> oraz współczynników niezależnych od rodzaju funkcji.

Wówczas







n

i

i

i

x

f

a

f

S

0

)

(

)

(

(1)





b

a

dx

x

f

f

I

)

(

)

(

(2)

Można wykazać, że suma ciągu S(f) jest zbieżna do całki oznaczonej I(f) dla





n

46

Uwagi ogólne o całkowaniu numerycznym

Całkowanie numeryczne stosujemy wtedy gdy trudno obliczyć całkę w sposób analityczny
lub w przypadku gdy znane są tylko wartości funkcji podcałkowej w punktach zwanych
węzłami

W celu przybliżonego obliczenia całki oznaczonej funkcji f(x) na skończonym przedziale

]

;

[ b

a

x



, funkcję podcałkową f(x) zastępujemy interpolującą funkcją φ(x), którą łatwo

można całkować. Funkcja φ(x) będzie wielomianem interpolującym z węzłami interpolacji x

0

,

x1, …, x

n

. Wówczas całkę możemy zastąpić sumą, współczynniki a

i

zależą od metody

obliczeń













b

a

b

a

n

i

i

i

x

f

a

dx

x

dx

x

f

0

)

(

)

(

)

(



(4)

Kwadratury z ustalonymi węzłami

Kwadratury Newtona-Cotesa
Całkowanie z zastosowaniem kwadratur o ustalonych węzłach polega na zastąpieniu funkcji
podcałkowej wielomianem interpolacyjnym Lagrange’a Li(x).
Jeżeli dla skończonego przedziału [a,b] wybierzemy zbiór punktów węzłowych {x

0

, ..., x

n

}

takich, że:

x

i

= x

0

+ i · h

47

gdzie:

n

a

b

h

)

(





dla wówczas=(0, 1, …, n)

n - oznacza również stopień wielomianu interpolacyjnego







n

i

i

i

x

L

x

f

x

0

)

(

)

(

)

(



(5), gdzie:













n

k

i

k

k

i

k

i

x

x

x

x

x

L

0

)

(

)

(

)

(

(6)

wówczas

 

   









 









 

7

dx

!

1

n

x

f

x

x

dx

x

L

x

f

dx

x

f

b

a

n

0

i

b

a

n

0

i

1

n

i

i

i

b

a























 













 





 

8

dx

x

f

x

x

!

1

n

1

dx

x

f

a

b

a

n

0

i

b

a

n

0

i

1

n

i

i

i





















ξ – dowolny punkt pośredni z przedziału (x,x

i

)

gdzie:

]

;

[

)

(

b

a

x





oraz





















b

a

n

i

n

i

i

i

t

n

t

t

t

i

n

i

h

dx

x

L

a

0

)

)...(

1

(

)!

(

!

)

1

(

)

(

(9)

ostatecznie wzór kwadratury Newtona – Cotesa oraz reszta kwadratury











b

a

n

i

i

i

f

E

x

f

a

dx

x

f

0

)

(

)

(

)

(

(10)













b

a

n

i

n

i

dx

x

f

x

x

n

f

E

0

)

1

(

))

(

(

)

(

)!

1

(

1

)

(



(11)

Twierdzenie 1

Kwadratury Newtona-Cotesa oparte na n+1 węzłach są rzędu











ych

nieparzyst

n

dla

1

n

parzystych

n

dla

2

n

wówczas wzory kwadratur dla n parzystych mają postać

48

 

 



 



 





  



 

 

12

b

,

a

gdzie

,

dt

n

t

...

1

t

t

!

2

n

f

h

dx

x

f

a

dx

x

f

b

a

n

0

i

n

0

2

n

2

n

3

n

i

i

b

a





























oraz dla nieparzystych mają postać:























b

a

n

n

i

n

n

n

i

i

dt

n

t

t

t

n

f

h

x

f

a

dx

x

f

0

0

)

1

(

)

2

(

)

)...(

1

(

)!

1

(

)

(

)

(

)

(



(13)

Kwadratury z ustalonymi węzłami – wzór trapezów

Wyznaczmy wzór kwadratury dla n = 1

 



   





  

 

 

14

f

E

dx

x

f

x

x

x

x

x

f

x

x

x

x

dx

x

f

1

0

x

x

1

0

1

0

0

1

0

1

b

a































gdzie:



















1

1

3

1

1

)

)(

(

)

(

2

1

)

)(

))(

(

(

2

1

)

(

x

x

x

x

o

n

o

n

o

o

dx

x

x

x

x

f

dx

x

x

x

x

x

f

f

E





)

(

12

2

)

(

3

)

(

2

1

3

0

1

1

0

2

0

1

3





n

n

f

h

x

x

x

x

x

x

x

x

x

f

























wówczas:































b

a

n

f

h

x

x

x

f

x

x

x

x

x

f

x

x

x

x

dx

x

f

)

(

12

)

(

)

(

2

)

(

)

(

)

(

2

)

(

)

(

3

0

1

1

0

1

2

0

0

1

0

2

1



49





)

(

12

)

(

)

(

2

)

(

3

1

0

0

1



n

f

h

x

f

x

f

x

x









(16)

Ostatecznie otrzymujemy

 

   





 

 

17

b

a

,

''

f

12

h

b

f

a

f

2

h

dx

x

f

3

b

a

















powyższe równanie jest znane jako wzór trapezów – widać, że wzór pozwala dokładnie
obliczyć całkę dla dowolnej funkcji, dla której druga pochodna jest zero (wtedy reszta jest 0).

Zauważmy, że dla powyższych rozwiązań końce przedziału całkowania są węzłami
kwadratury x

0

=a oraz x

n

=b. Wzory kwadratur oparte na tych węzłach nazywamy wzorami

zamkniętymi Newtona – Cotesa.
Podobnie wyznaczyć można wzory dla innych n.

Kwadratury Newtona-Cotesa – wzory zamknięte

Poniżej znajdują się wzory zamknięte dla n = (1..3):

wzór trapezów

 

   





 

 

18

,

''

12

2

1

1

0

3

1

0

x

x

f

h

x

f

x

f

h

dx

x

f

n

b

a



















wzór parabol (Simpsona)

 

 

   





 

 

19

,

90

4

3

2

2

0

)

3

(

5

2

1

0

x

x

f

h

x

f

x

f

x

f

h

dx

x

f

n

b

a





















wzór Bessela

 

 

 

   





 

 

20

,

80

3

3

3

8

3

3

3

0

)

4

(

7

3

2

1

0

x

x

f

h

x

f

x

f

x

f

x

f

h

dx

x

f

n

b

a























50

Interpretacja graficzna

wzór parabol (Simpsona)

wzór Bessela

Kwadratury Newtona - Cotesa - wzory otwarte

Wzory kwadratur nie opartych na węzłach będących końcami przedziału całkowania
nazywamy wzorami otwartymi Newtona - Cotesa

51

Jeżeli dla skończonego [a ; b] przedziału wybierzemy zbiór punktów węzłowych {x

0

, …, x

n

}

takich, że :

h

i

x

x

i







0

gdzie:

2

)

(







n

a

b

h

dla i=(0, …, n)

oraz

h

a

x





0

W oparciu o Twierdzenie 1 można również wyznaczyć wzory kwadratur otwartych.

Wówczas

wzory kwadratur dla n parzystych mają postać

 

 





 





  



 

21

]

b

,

a

[

gdzie

,

dt

n

t

...

1

t

t

!

2

n

f

h

dx

x

f

a

dx

x

f

n

0

i

1

n

1

2

2

n

3

n

i

i

b

a

































oraz dla n nieparzystych mają postać:































b

a

n

i

n

n

n

n

i

i

b

a

gdzie

dt

n

t

t

t

n

f

h

x

f

a

dx

x

f

0

1

1

)

1

(

)

1

(

)

2

(

)

22

(

]

;

[

,

)

)...(

1

(

)!

1

(

)

(

)

(

)

(





wzór prostokątów

 

 

 

 

23

x

x

,

''

f

3

h

x

f

2

h

dx

x

f

0

n

1

1

3

0

b

a



















 

   





 

 

24

x

x

,

''

f

3

h

3

x

f

x

f

2

h

3

dx

x

f

1

n

2

1

3

1

0

b

a



















 

   

 





 

 

25

,

45

14

2

2

3

4

2

3

1

)

4

(

5

2

1

0

x

x

f

h

x

f

x

f

x

f

h

dx

x

f

n

b

a























 

     

 





 

 

26

,

144

95

11

11

24

5

3

4

1

)

4

(

5

3

2

1

0

x

x

f

h

x

f

x

f

x

f

x

f

h

dx

x

f

n

b

a

























52

Interpretacja geometryczna

wzór prostokątów

wzór trapezów

Kwadratury Newtona-Cotesa – porównanie

Przykładowe wyniki obliczeń całki

 

 





29289922

,

0

2

2

1

x

cos

dx

x

sin

4

0

4

0

















Wyniki obliczeń dla wzorów otwartych i zamkniętych:

n

0

1

2

3

wzory zamknięte

0,27768018

0,29293264

0,29291070

błąd

0,01521303

0,00003924

0,00001748

wzory otwarte

0,30055887

0,29798754

0,29285866

0,29286923

błąd

0,007666565

0,00509432

0,00003456

0,00002399

53

Kwadratury Newtona-Cotesa – podsumowanie

Kwadratury Newtona-Cotesa dają coraz lepszą dokładność wraz ze wzrostem n. Jednak wraz
ze wzrostem n wzór na sumę posiada również rosnącą liczbę czynników. Dlatego kwadratur
Newtona-Cotesa nie stosuje się dla dużych n oraz z uwagi na pojawianie się ujemnych
współczynników (dla n = 7). Co więcej dla dużych n istnieją funkcje ciągłe, dla których ciąg
kwadratur Newtona – Cotesa nie jest zbieżny do całki



b

a

dx

x

f

)

(

W praktyce nie stosuje się kwadratur Newtona – Cotesa wysokiego rzędu.
Jak widać aby obliczyć wartość całki stosując kwadratury Newtona – Cotesa wystarczy
wyznaczyć wartość funkcji w punktach węzłowych, a następnie podstawić je do wzoru.
Większą trudność sprawia oszacowanie reszty kwadratury, gdyż należy wyznaczyć
maksimum pochodnej k rzędu funkcji dla przedziału całkowania.

Obliczyć poniższą całkę stosując metodę Simpsona

 

 





4

4

5

0

d

)

sin(

d

2

90

1

sin

2

sin

4

0

sin

6

dx

)

x

sin(















 





























 













 











cos

d

)

sin(

d

4

4

, wtedy





1

)

cos(

max













,

0

wówczas wyniki:

09439510

,

2

3

2

dx

)

x

sin(

0











oraz oszacowanie błędu obliczeń:

106256835

,

0

2880

1

2

90

1

5

2























 



rozwiązanie dokładne:





2

)

0

cos(

)

cos(

)

x

cos(

dx

)

x

sin(

0

0





















Kody źródłowe:

przykładowa całkowana funkcja:

function f(x:extender):extender,
begin
result:=sin (2*x*x)+2*cos(3x);
end;

54

wzór prostokątów:

function calka_n_0_wzor_otwarty(a,b:extended):extended;
var h:extended;
begin
h:=(b-a)/2;
result:=2*h*f(a+h);
end;


wzór trapezów:

function calka_n_1_wzor_zamkniety(a,b:extended):extended;
var h:extended;
begin
h:=(b-a);
result:=(h/2)*(f(a)+f(b));
end;

wzór parabol:

function calka_n_2_wzor_zamkniety(a,b:extended):extended;
var h:extended;
begin
h:=(b-a)/2;
result:=(h/3)*(f(a)+f(a+h)=f(b));
end;

55

Kwadratury z ustalonymi węzłami

Złożone kwadratury Newtona – Cotesa:
W poprzednim punkcie zostały przedstawione kwadratury Newtona-Cotesa. Wiemy, że błąd
kwadratury zależy od długości przedziału całkowania oraz rzędu kwadratury. Dlatego w
oparciu właśnie o kwadratury Newtona-Cotesa niskich rzędów zdefiniujemy tzw. wzory
złożone Newtona-Cotesa. W celu zwiększenia dokładności obliczeń należy zmniejszyć
długość przedziału. Wynika z tego, że dzieląc przedział całkowania na m podprzedziałów i
obliczając dla każdego przedziału kwadraturę Newtona-Cotesa niskiego rzędu, następnie, gdy
zsumujemy wyniki podprzedziałów, uzyskamy wynik zbieżny do rozwiązania dokładnego,
natychmiast zwiększy się dokładność obliczeń.

m=1

Otrzymane w ten sposób kwadratury określono na całym przedziale całkowania nazywamy
złożonymi kwadraturami Newtona-Cotesa

Twierdzenie 2:
Podzielmy przedział całkowania [a,b] na m części przy pomocy punktów:

h

i

a

x

i







, dla

)

1

m

,...,

1

,

0

(

i





,

m

)

a

b

(

h





mamy wtedy:

 











1

m

0

i

x

x

b

a

1

i

i

dx

)

x

(

f

dx

)

x

(

f

56

stosując do każdej z całek po prawej stronie równania dowolny wzór kwadratury Newtona-
Cotesa

...

Na podstawie poprzedniej zależności możemy wyprowadzić wzory złożonych kwadratur:

złożony wzór prostokątów:

E

2

h

x

f

h

E

2

h

x

f

h

dx

)

x

(

f

1

m

0

i

i

1

m

0

i

i

1

m

0

i

i

b

a

























































(28)

gdzie błąd przybliżenia

 











1

m

0

i

)

2

(

3

f

3

h

E

(29)

oraz





1

i

i

i

x

,

x







złożony wzór trapezów:





 

E

x

f

2

)

x

(

f

)

x

(

f

h

E

)

x

(

f

)

x

(

f

2

h

dx

)

x

(

f

1

m

1

i

i

m

0

1

m

0

i

i

1

m

0

i

1

i

i

b

a

















































(30)

gdzie błąd przybliżenia

 

i

1

m

0

i

)

2

(

3

f

12

h

E













(31)

oraz





1

i

i

i

x

,

x







57

złożony wzór parabol – dla m parzystych:









































1

m

0

i

i

1

m

0

i

1

i

1

i

i

b

a

E

)

x

(

f

)

2

h

x

(

f

4

)

x

(

f

3

h

dx

)

x

(

f

(32)

E

)

x

(

f

4

)

x

(

f

2

)

x

(

f

)

x

(

f

3

h

dx

)

x

(

f

1

2

m

1

i

1

i

2

1

2

m

1

i

i

2

m

0

b

a















































(33)

gdzie błąd przybliżenia

)

(

f

h

180

)

a

b

(

E

)

4

(

4









(34)

oraz





1

i

i

i

x

,

x







Podsumowanie:

Z powyższych wzorów wynika, że jeśli funkcja f jest wystarczająco regularna, to całka tej
funkcji może być z dowolnie dużą dokładnością przybliżona przy pomocy złożonych
kwadratur. Dla złożonych kwadratur Newtona-Cotesa, również największą trudność sprawia
oszacowanie błędu obliczeń.

Twierdzenie 3:
Jeżeli

]

b

;

a

[

c

f



(jest jednokrotnie całkowalna w przedziale [a,b]), to dla każdego ciągu m

złożonych kwadratur Newtona-Cotesa n-tego rzędu dla funkcji f jest zbieżny do całki:



b

a

dx

)

x

(

f

, przy





m

inaczej

 















1

m

1

i

x

x

m

b

a

1

i

i

dx

)

x

(

f

lim

dx

)

x

(

f

h

i

a

x

i







,

)

1

m

,...,

1

,

0

(

i





,

m

)

a

b

(

h





58

Przykład:
Obliczyć poniższą całkę:

Rozwiązanie dokładne





59815

,

53

0

4

4

0









e

e

dx

e

x

Metoda Simpsona m=1, h=2





76958

,

56

4

3

2

4

2

0

4

0











e

e

e

dx

e

x

błąd = -3,17143

Metoda Simpsona m=1, h=1



 



86385

,

53

4

3

1

4

3

1

4

3

2

2

1

0

4

2

2

0

4

0

























e

e

e

e

e

e

dx

e

dx

e

dx

e

x

x

x

błąd = -0,26570


Metoda Simpsona m=4, h=0,5

61622

,

53

4

3

3

2

2

1

1

0

4

0

























dx

e

dx

e

dx

e

dx

e

dx

e

x

x

x

x

x

błąd = -0,01807


Przykładowe wyniki obliczeń całki:





2

)

cos(

)

sin(

0

0













x

dx

x

Wyniki obliczeń dla wzorów złożonych w zależności od m:

m

2

4

8

16

32

64

wzór

prostokątów

wynik

1,110720735 1,751785441 1,936297295 1,983970758 1,995986709 1,998996147

błąd

0,889279265 0,248214559 0,066102705 0,016029242 0,004013791 0,001003853

wzór

trapezów

wynik

1,570796327 1,896118898 1,1974231602 1,993570344 1,998393361 1,999598389

błąd

0,429203673 0,103881102 0,025768398 0,006429656 0,001606639 0,000401611

wzór parabol

wynik

2,094395102 2,004559755 2,000269170 2,000016591 2,000001033 2,000000065

błąd

-0,094395102 -0,004559755 -0,000269170 -0,000016591 0,000001033 0,000000065

Kody źródłowe:

Złożony wzór prostokątów:

function
Calka_Prostokat_Zlozony(a,b:extended;m:integer):extended;
var h,x,I_0:extended;
begin
h:=(b-a)/m;

59

I_0:=0.0;
for i:=0 to m-1 do
begin
x:=a+(i*h);
x:=x+(h/2);
I_0:=I_0+f(x)
end;
result:=I_0*h
end;

Złożony wzór trapezów:

function
Calka_Trapez_Zlozony(a,b:extended;m:integer):extended;
var h,x,I_0,I_1:extended;
begin
h:=(b-a)/m;
I_0:=f(a)+f(b);
I_1:=0.0;
for i:=0 to m-1 do
begin
x:=a+(i*h);
I_1:=I_1+f(x);
end;
result:=(0.5*I_0+I_1)*h;
end;


Złożony wzór parabol:

function
Calka_Simpson_Zlozony(a,b:extended;m:integer):extended;
var h,x,I_0,I_1,I_2:extended;
begin
h:=(b-a)/m;
I_0:=f(a)+f(b);
I_1:=0.0;
I_2:=0.0;
for i:=0 to m-1 do
begin
x:=a+(i*h);
if (i mod 2)=0
then I_2:=I_2+f(x);
else I_1:=I_1+f(x);
end;
result:=(I_0+2.0*I_2+4.0*I_1)*h/3;
end;

60

Kwadratury z ustalonymi węzłami

Metoda ekstrapolacyjna – Romberga

Metoda adaptacyjna

Kwadratury Gaussa

W poprzednim rozdziale omawiane były metody Newtona-Cotesa, które bazowały na
całkowaniu wielomianu interpolacyjnego stopnia n. Pozwalał on na dokładne przybliżenie
całki dla funkcji stopnia (n+1) lub (n+2). Wszystkie metody Newtona-Cotesa
wykorzystywały wartości funkcji dla punktów węzłów równoległych. Ten warunek jest
dogodny dla konstrukcji złożonych wzorów, ale może znacząco zmniejszyć dokładność
przybliżenia całek.
Poniżej widać kwadraturę Newtona-Cotesa opartą na końcach przedziału – przybliżenie
funkcji nie jest dokładne.

Kwadratury Gaussa dobierają takie węzły, aby osiągnąć optymalne przybliżenie całki:









n

1

i

i

i

b

a

)

x

(

f

c

dx

)

x

(

f

Węzły kwadratury x

1

, x

2

,..., x

n

z przedziału całkowania [a,b] oraz współczynniki c

1

, c

2

,...,c

n

są

tak dobrane, aby zminimalizować błąd przybliżenia. Nie zakładamy jednak żadnych
ograniczeń na węzły x

i

o współczynnikach c

i

natomiast chcemy zmaksymalizować rząd

kwadratury.

61

Mamy znaleźć 2n stałe we wzorze, przypuszczamy więc, że rząd kwadratury powinien
wynosić 2n, a więc wzór powinien być dokładny dla wielomianów stopnia (2n-1). Dla
wielomianu stopnia równego lub wyższego od (2n-1) należy uwzględnić resztę kwadratury E:

E

)

x

(

f

c

dx

)

x

(

f

n

1

i

i

i

b

a











Aby pokazać proces wyboru odpowiednich węzłów oraz współczynników zilustrujemy to na
przykładzie:
n=2 – liczba węzłów oraz przedział [-1,1]

Należy wyznaczyć x

1

, x

2

oraz c

1

, c

2

, ponieważ:

)

x

(

f

c

)

x

(

f

c

dx

)

x

(

f

2

2

1

1

1

1









(64)

Zgodnie z założeniami, powyższe przybliżenie będzie dokładne tylko gdy stopień funkcji f(x)
mniejszy lub równy

3

1

)

2

(

2

1

2











n

wówczas

3

3

2

2

1

0

x

a

x

a

x

a

a

)

x

(

f









(65)

gdzie a

0

, a

1

, a

2

, a

3

– dowolne stałe

a więc możemy zapisać:





























dx

x

a

dx

x

a

xdx

a

dx

1

a

dx

x

a

x

a

x

a

a

3

3

2

2

1

0

3

3

2

2

1

0

Poszukujemy x

1

, x

2

oraz c

1

, c

2

, więc można zapisać:

2

1

c

1

c

2

1









(67)

0

dx

x

x

c

x

c

1

1

2

2

1

1















(68)

3

2

dx

x

x

c

x

c

1

1

2

2

2

2

2

1

1















(69)

0

dx

x

c

x

c

1

1

3

2

2

3

1

1















(70)

otrzymaliśmy układ liniowych równań algebraicznych, którego rozwiązanie jest:

1

c

1



1

c

2



3

3

x

1





3

3

x

2



62

ostatecznie otrzymaliśmy wzór na przybliżenie całki dla funkcji danej wzorem (64):



































3

3

f

3

3

f

dx

)

x

(

f

1

1

Metoda poszukiwania węzłów oraz współczynników może zostać zastosowana dla funkcji
wyższych rzędów, lecz istnieje łatwiejsza metoda. Bazuje ona na wykorzystaniu jednego z
wielomianów ortogonalnych. Mają one taką właściwość, że całka z iloczynu dwóch
dowolnych wielomianów jest równa 0:













b

a

dx

x

g

x

f

g

f

0

)

(

)

(

,

Dla naszych dalszych rozważań, zajmiemy się tylko ciągiem wielomianów Legendre’a









)

(

)

(

),

(

1

0

x

P

x

P

x

P

n

, który posiada właściwości:

1. Dla każdego n, P

n

(x) jest wielomianem stopnia n

2.









1

1

0

)

(

)

(

dx

x

P

x

P

n

, gdzie P(x) jest wielomianem stopnia wyższego

(niższego???)

od n

Wzory kilku pierwszych wielomianów Legendre’a:

1

)

(

0



x

P

x

x

P



)

(

1

















3

1

2

3

)

(

2

2

x

x

P











 



x

x

x

P

3

5

2

5

)

(

3

3

oraz wzór rekurencyjny do wyznaczania wielomianów krzywych

???

oraz ich pochodnych:

)

(

1

)

(

1

2

)

(

2

1

x

P

n

n

x

P

x

n

n

x

P

n

n

n



















)

(

1

)

(

1

)

(

'

1

2

2

x

P

x

n

x

P

x

x

n

x

P

n

n

n

















Pierwiastki tego wielomianu są tylko pojedyncze, leżą w przedziale <-1,1> oraz są
symetryczne względem początku układu współrzędnych. Najważniejszą właściwością jest to,
że bardzo dobrze nadają się do wyznaczania punktów węzłowych oraz współczynników.

Punkty węzłowe potrzebne w celu dokładnego przybliżenia całki funkcji stopnia mniejszego
od 2n są pierwiastkami wielomianu Legendre’a

n-tego stopnia.

Jeżeli x

1

, x

2

, ... x

n

są pierwiastkami wielomianu Lagrange’a, P

n

(x) stopnia n dla każdego i =

(1,2,...n) wtedy współczynniki c

i

można obliczyć ze wzoru:









2

1

2

1

1

1

1

1

)

(

)

1

(

1

2

)

(

'

)

(

)

1

(

2

i

n

i

i

n

i

n

n

i

k

k

k

i

k

x

P

n

x

x

P

x

P

n

dx

x

x

x

x

c



































wówczas jeżeli P(x) jest dowolnym wielomianem stopnia mniejszego od 2n wtedy:













1

1

1

)

(

)

(

n

i

i

i

x

P

c

dx

x

P

(75)

63

reszta dla stopnia większego od (2n-1):









 













)

2

2

(

1

1

2

!

1

2

)

(

n

f

n

dx

x

P

E

Pierwiastki wielomianu oraz współczynniki c

i

są tabelaryzowane – nie obliczamy ich przy

każdym całkowaniu.

Prawa strona wzoru (75) jest nazwana kwadraturą Gaussa-Lagrange’a.

Wzór ten pozwala na obliczenie całki tylko w przedziale

1

,

1



. W celu obliczenia

przybliżenia całki w dowolnym przedziale

b

a,

dla którego f(x) jest ciągła należy wykonać

przekształcenie poprzez podstawienie.
Dokonujemy przekształcenia dowolnego przedziału

b

a,

do przedziału

1

,

1



.

2

2

a

b

t

a

b

x











dt

a

b

dx

2





Wówczas możemy zapisać:































b

a

dt

a

b

t

a

b

f

a

b

dx

x

f

1

1

2

2

2

)

(

po zastosowaniu kwadratury Gaussa-Lagrange’a otrzymujemy





























































2

2

2

2

2

2

)

(

1

1

1

a

b

x

a

b

f

c

a

b

dt

a

b

t

a

b

f

a

b

dx

x

f

i

n

i

i

b

a

gdzie x

i

– pierwiastek wielomianu stopnia n

c

i

– współczynnik ze wzoru (74)

przykład:
Obliczyć całkę wykorzystując kwadraturę Gaussa:







5

,

1

1

1093643

,

0

2

dx

e

x

dokończyć ???

najpierw przekształcenia przedziału całkowania:

dt

t

f

dt

t

f

dx

e

x



































































1

1

5

,

1

1

1

1

2

16

5

4

1

2

1

5

,

1

2

1

5

,

1

2

1

5

,

1

2

po zastosowaniu kwadratury Gaussa-Lagrange’a, n=2 otrzymujemy:

c

1

=1 c

2

=1

3

3

1





x

3

3

2



x

64

1094003

,

0

1

1

4

1

4

1

16

5

3

3

16

5

3

3

16

5

2

16

5

1

5

,

1

1

2

1

2































































































































































































e

e

e

c

e

c

dx

e

x

x

x

Po zastosowaniu kwadratury Gaussa-Lagrange’a, n=3 otrzymujemy:

9

5

1



c

9

8

2



c

9

5

3



c

3

15

1



x

0

2



x

3

15

3





x

 

1093642

,

0

9

5

9

8

9

5

4

1

4

1

16

5

3

15

16

5

0

16

5

3

15

16

5

3

16

5

2

16

5

1

5

,

1

1

3

2

1

2

























































































































































































































e

e

e

e

c

e

c

e

c

dx

e

x

x

x

x

Tablica pierwiastków wielomianu Legendre’a:

n

x

0

x

1

x

2

x

3

x

4

1

0

-

-

-

-

2

-0,5773509

0,5773509

-

-

-

3

0,7745966

0

-0,77459666

-

-

4

0,86113631

0,33998104

-0,33998104

-0,86113631

-

5

-0,53846931

-0,90617984

0

0,90617984

0,53846931

Kod źródłowy kwadratury Gaussa:
Przykład f. całką wykorzystując kwadraturę G w przedziale

b

a,

oraz n- oznacz. stopnia

met. L.


function Gauss_Legendre (n:Integer; var x:vector); Extended
var k,l,m:Integer;
p0,p1,p2,s,x1,ci:Extended
begin
if n>0 then
begin
s:=0
l:=n div2;
for k:=0 to n do
begin
p0:=1;
x1:=x[k];

65

p1:=x1
for m:=1 to n-1 do
begin
p2:=((2*m+1)*x1*p1-m*p0)/(m+1);
p0:=p1;
p1:=p2;
end;
ci:=-2*(x1*x1-1)/((n+1)*(n+1)*p1*p1);
s:=s+f(x1)*ci;
end;
Result:=s;
end;
end;

gdzie x – tabl. zawiera pierwiastki wielomianu Legendre’a
oraz n – stopień wielomianu Legendre’a

Obliczanie całek niewłaściwych:
Przy numerycznym całkowaniu w przedziale obustronnie lub jednostronnie nieskończonym,
należy stosować takie funkcje wagowe, które zapewniają zbieżność całki, W zależności od
przedziału stosujemy odpowiedni wielomian interpolacyjny:













n

i

i

i

b

a

x

P

c

dx

x

f

x

w

I

1

)

(

)

(

)

(

Całkowanie funkcji osobliwych:
W celu obliczenia przybliżonej całki funkcji posiadającej punkty nieciągłości, należy
zmodyfikować poznanie metody. Aby dokonać takiego całkowania, należy zlokalizować
wszystkie punkty nieciągłości należące do przedziału całkowania.

Przykłady funkcji nieciągłych:





b

a

dx

x

f

I

)

(

,

gdzie f(x):

przypadek 1*

przypadek 2*

0





66











b

a

c

a

b

a

dx

x

f

c

dx

x

f

dx

x

f

)

(

)

(

)

(

Powyższe metody wymagają wielu przekształceń matematycznych, które nie dają wyznaczyć
się informatycznie.
Ma to charakter historyczny.

Alternatywną metodą wyliczenia całki jest metoda polegająca na oddalaniu się od punktu
nieciągłości

o małe ε ???

Metoda zwykłego szeregu Taylora:

9253141

,

2

0017691

,

0

9235450

,

2

1

0









dx

x

e

x

Metoda adaptacyjna:

925324091

,

2

1

0000000001

,

0

1

0









dx

x

e

dx

x

e

x

x

zał. dokładne 0,0001 liczba przedzi. 235 il. wywoł. funkcji 941

???

Całki wielokrotne:
rozpatrzmy całkę podwójną:

 

b

a

x

d

x

c

dx

dy

y

x

f

)

(

)

(

)

,

(

x – stała i stosunek kwadratury do obliczonej całki





)

(

)

(

)

,

(

)

(

x

d

x

c

dy

y

x

f

x

s

67

otrzymujemy:



































1

1

1

1

2

2

2

0

)

,

(

4

)

,

(

2

)

,

(

)

,

(

6

)

(

)

(

)

(

m

j

m

j

j

j

m

y

x

f

y

x

f

y

x

f

y

x

f

m

x

c

x

d

x

s

m – liczba przedziałów przedziału [c(x) ; d(x)]

wykorzystanie z kolei kwadratury Simpsona do obliczenia całki:





b

a

dx

x

S

s

)

(

ostatecznie:



































1

1

1

1

2

2

2

0

)

(

4

)

(

2

)

(

)

(

6

n

i

n

i

i

i

n

x

S

x

S

x

S

x

S

n

a

b

s

np.:









5

,

1

0

,

1

0

,

2

4

,

1

4295549269

,

0

)

2

ln(

dx

dy

y

x

stos. złoż. kwadratury Simpsona:
dla n=4 i m=2









5

,

1

0

,

1

0

,

2

4

,

1

7

4295522438

,

0

)

2

ln(

dx

dy

y

x



 





1800

6

,

0

5

,

0





E

6

10

72

,

4









E

Całkowanie wielokrotne – całka podwójna
kod źródłowy:
function f (x,y:real):real;
begin
Result:=

explg (x)???

end;

function c(x:real):real;
begin
Result:=x*x*x
end;

function d(x:real):real;
begin
Result:=x*x
end;

68

function simpson-Podwojna
(a,b:real; m,n:integer):real;
var i,j:integer
h,hx,j1,j2,j3
k1,k2,k3,x,y,x:real;
begin
h:=(b-a)/(2*n)
j1:=0; j2:=0;
j3:=0
for i:=0 to 2*n do
begin
x:=a+i*h;
nx:=(d(x)-c(x))/(2*m);
k1:=f(x,c(x))+f(x,d(x));
k2:=0
k3:=0
for j:=1 to 2*m-1 do
begin
y:=c(x)+j*hx
z:=f(x,y);

if Odd(j) ???

then k3 :=k3+z
else k2:=k2+z
end
i:=(k1+2*k2+4*k3)+hx/3
if(i=0) or (i=2*n)
then j1:=j1+1

else if Odd (i)
then j3:=j3+l
else j2:=j2+l
end;
Result:=(j1+2*j2+4*h/3)
end;

Kody źródłowe:

Function f(x,y:real):real;
Begin
Result:=exp(y/x);
End;

Function c (x:real):real;
Begin
Result:=x*x*x;
End;

Function d (x:real):real;
Begin
Result:=x*x;

69

End;

function Simpson_Podwojna (a,b:real; m,n:integer):real;
var i,j:integer;
h,hx,j1,j2,j3,k1,k2,k3,I,x,y,z:real;
begin
h:=(b-a)/(2*n);
j1:=0;
j2:=0;
j3:=0;
for i:=0 to 2*n do
begin
x:=a+i*h;

tu tomek napisał ... j*h a ja mam i*h sprawdźcie!!!

hx:=(d(x)-c(x))/(2*m)
k1:=f(x,c(x)+f(x,d(x));
k2:=0;
k3:=0;
for j:=1 to 2*m-1 do
begin
y:=c(x) +j*hx;
z:=f(x,y);
if odd(j)
then k3:=k3+z
else k2:=k2+z
end;
I:=(k1+2*k2+4*k3)*hx/3;
If (i=0) or (i=2*n)
then j1:=j1+I

else if Odd(i)

tu też mam troche inaczej

then j3:=j3+l
else j2:=j2+l
end;
Result:=(j1+2*j2+4*j3)*h/3;
End;

Obliczanie całek metodą Monte-Carlo

Metoda Monte Carlo – całka pojedyncza

Mamy obliczyć przybliżoną wartość In całki oznaczonej





a

b

f(x)dx

I

(129)

przy założenie, że f(x) jest funkcją ciągłą w przedziale domkniętym [a;b]
Następnie n-krotnie generujemy realizację x zmiennej losowej X o rozkładzie równomiernym
w przedziale [a,b], otrzymujemy w rezultacie ciąg liczb x

1

, x

2

, ..., x

n

,

70

Przybliżoną wartość całki określa wzór









n

i

i

n

x

f

n

a

b

I

1

)

(

)

(

(130)

Zbieżność przybliżonej wartości całki In do jej dokładnej wartości jest gwarantowana w
sensie probalistycznym, tzn., że dla każdego ε>0 zachodzi relacja:

1

)

1

n

(

P

lim

n















(131)

gdzie P jest prawdopodobieństwem, przy czym błąd metody monte Carlo można określić
wzorem

n

I

I

n

1





(132)

Maszyny typowe zazwyczaj dysponują generatorem liczb losowych o rozkładzie
równomiernym w przedziale [0,1], zachodzi zatem konieczność przeliczania wylosowanej
realizacji Y z przedziału [0,1], na realizację zmiennej losowej X o rozkładzie równomiernym
w przedziale [a, b]. Wzory poniżej, to ilustrują:
Przeliczanie zmiennych losowych Y na X

X=(b-a)Y+a

(133)

Przeliczanie realizacji zmiennych losowych Y na X

x=(b-a)y+a

(134)

Kod źródłowy:

Function Calka_MonteCarlo (a,b:extended; n:lingint):extended
Var i:integer;
h,x:extended;

begin
Result:=0;
h:=(b-a)/n;
for :=1 to n do
begin
x:=a+((b-a)*Random);
Result:=Result+f(x);
End;
Result:=Result*h;
end;

71

Całkowanie – Metody Monte Carlo. Całki wielokrotne właściwe

Dana jest całka wielokrotna:

3

2

1

2

1

...

)

,...,

,

(

dx

dx

dx

x

x

x

f

I

n







(135)

Każda ze zmiennych x

i

zawiera się w przedziale

i

i

i

b

x

a





(136)

Losujemy n-krotnie punkty Pk kostki, przyjmując że współrzędnymi tych punktów są
zmienne losowe X

1

(k)

, X

2

(k)

,..., X

m

(k)

o rozkładzie równomiernym w odpowiednich

przedziałach [a

i

,b

i

]. Każdemu punktowi Pk odpowiada zatem ciąg m-elementowy

wylosowanych realizacji x

1

(k)

, x

2

(k)

,..., x

m

(k)

. Po wylosowaniu sprawdzamy czy każdy punkt Pk

należy do obszaru Ω. Możemy zatem po zakończeniu losowań uporządkować wylosowane
punkty w następujący sposób:





k

P

dla k=(1, 2, ..., n),





k

P

dla k=(n+1, n+2, ..., N)

(137)

Przy dostatecznej liczbie losowań możemy obliczyć przybliżoną wartość V

Ω

objętości m-

wymiarowanego obszaru całkowania:

V

N

n

V





(138),

gdzie

)

(

1

i

m

i

i

a

b

V









(139)

jest objętością m-wymiarowej kostki. Przybliżona wartość I

n

całki (135) określana wzorem:





n

k

k

n

P

f

N

V

I

1

)

(

(140)

Aby uzyskać ciąg realizacji X1, X2, ..., Xm losujemy m-krotnie zmienną losową Y o
rozkładzie równomiernym w przedziale [0,1], a następnie przeliczamy realizacje y na
realizację x

i

(k)

w przedziale [a

i

; b

i

] korzystając ze zależności:

x

i

(k)

=(b

i

-a

i

)y

i

(k)

+a

i

Rozwiązywanie równań różniczkowych
Rozważania dotyczą:

- równań różniczkowych zwyczajnych z warunkami początkowymi i brzegowymi,
- układów równań różniczkowych zwyczajnych z warunkami początkowymi i

brzegowymi,

- równań różniczkowych cząstkowych,
- układów równań różniczkowych cząstkowych

72

Numeryczne obliczanie pochodnej
Aby obliczyć pochodną funkcji f w punkcie x

0

, mamy

h

x

f

h

x

f

x

f

h

)

(

)

(

lim

)

(

'

0

0

0











(1)

W celu wyznaczenia przybliżonej wartości pochodnej f’(x) dla małych wartości h możemy
użyć wyrażenia

h

x

f

h

x

f

)

(

)

(

0

0





(2)

Możemy to udowodnić. Dokonujemy interpolacji f(x); zdefiniujemy wielomian Lagrange`a
stopnia 1, który będzie oparty na węzłach x

0

oraz x

0

+h pod warunkiem, że funkcja f jest

dwukrotnie różniczkowalna w przedziale [a; b]

!

2

)

)(

(

)

(

)

(

0

0

1

h

x

x

x

x

x

P

x

f











ktoś nie zdązył do końca, a ja nie wiem czy w

mianowniku jest 2!, czy co??

Wówczas obliczając pochodną funkcji f, mamy:











h

)

x

(

f

)

h

x

(

f

)

x

(

'

f

0

0

0

(3)

oraz po uproszczeniu

)

(

"

2





f

h





(4)

73

Jeżeli dokonamy interpolacji funkcji f wielomianem stopnia 2, wówczas













h

h

x

f

h

x

f

x

f

2

)

(

)

(

)

(

'

0

0

0

(5)

oraz po uproszczeniu

)

(

6

)

3

(

2





f

h





I(6)

Za pomocą wielomianów wyższych stopni możemy wyznaczyć wzory na przybliżenie
pochodnej.

Przykład: Obliczyć pochodną funkcji f(x)=ln(x) dla x

0

=1,8 stosując metodę w przód i wstecz

stopnia 1.

2

)

(

)

(

)

(

'

0

0

0

h

h

x

f

h

x

f

x

f









M

h

h

f

h

f

f

2

)

8

,

1

(

)

8

,

1

(

)

(

"











,

gdzie:

)

(

"

max

)

8

,

1

,

8

,

1

(





f

M

h







brak textu

Metoda w przód

Metoda wstecz

h

Wynik

Błąd

Wynik

Błąd

0,1

0,5406772

0,0154321

0,5715841

0,0173010

0,01

0,5540180

0,0015432

0,5571045

0,0015605

0,001

0,5554013

0,0001543

0,5557099

0,0001545

74

Numeryczne obliczanie drugiej pochodnej.
Jeżeli dokonamy rozwinięcia funkcji f w szereg Taylora 3 rzędu w otoczeniu punktu x

0

dla

punktów x

0

+h oraz x

0

-h, zakładając, że f jest 4-krotnie różniczkowalna na przedziale [x

0

-h,

x

0

+h], wówczas dla x

0

-h:

Wzor (8)

gdzie

]

x

h;

[x

ξ

0

0





dodając oba równania stronami otrzymujemy:

wzór (9)

Rozwiązując równanie ze względu na f”(x

0

), otrzymujemy:















)]

(

)

(

2

)

(

[

1

)

(

"

0

0

0

2

0

h

x

f

x

f

h

x

f

h

x

f

(11)

oraz

)

(

[

12

)

(





n

n

f

h



, gdzie

]

;

[

0

0

h

x

h

x









Powyższy wzór otrzymamy różniczkując obustronnie wzór na pierwszą pochodną (5) kroku
h/2.

Przykład: Obliczyć drugą pochodną funkcji f(x)=ln(x

3

) dla x

0

=2 stosując powyższą metodę

różnych wartości kroku h.

Wzór

Rozwiązanie dokładne

2

3

)

(

"

x

x

f





f”(2)=-0,1875

h

Wynik

Błąd

Oszacowanie ε

0,1

-0,1875586182

0,0000586182

0,01

-0,1875005859

0,001

-0,1875000049

Równania różniczkowe zwyczajne z warunkami początkowymi

Równania różniczkowe są wykorzystywane do odwzorowania wielu problemów z zakresu
różnych nauk, które dotyczą zmian jednej zmiennej względem innej zmiennej. Wiele
problemów świata rzeczywistego opisywanych przez równania różniczkowe jest za bardzo
złożonych, aby uzyskać rozwiązanie dokładne (analityczne), należy znaleźć rozwiązanie
przybliżone. Są dwa sposoby aby uzyskać rozwiązanie przybliżone. Po pierwsze można
dokonać uproszczenia zadanego równania różniczkowego do takiej postaci, dla której znamy
rozwiązanie dokładne, wówczas uproszczone równanie będzie przybliżało rozwiązanie
równania zadanego. Drugim sposobem jest zastosowanie jednej z metod numerycznych w
celu przybliżenia rozwiązania zadanego równania.

Metody różnicowe jednokrokowe:

- metody Eulera;
- metody Rundego-Kutty
- metody Rundego-Kutty-Fehenberga

75

Metody różnicowe wielokrokowe:

- metody bezpośrednie
- metody pośrednie
- metoda Predictor-Corrector
- metoda trapezów

Metoda Eulera
Rozpatrzmy zadanie znalezienia funkcji y=y(t), które dla

]

,

[

b

a

t



spełnia równanie

różniczkowe zwyczajne pierwszego rzędu:

)

,

( y

t

f

dx

dy



(7)

oraz warunek początkowy

y(a)=y

0

(8),

gdzie f jest daną funkcją dwu zmiennych. Warunek wystarczający istnienia i jednoznaczności
rozwiązania zadania (7) oraz (8) jest warunek ciągłości funkcji f(t,y)

Metoda ta rzadko jest używana w praktyce, natomiast ze względu na jej prostotę, zostanie
przedstawiona w celu zobrazowania i zrozumienia technik stosowanych w bardziej
zaawansowanych metodach, pomijając przy tym złożone działania matematyczne.
Przedmiotem rozważań będzie poniższe równanie.

)

,

( y

t

f

dx

dy



, gdzie

]

,

[

b

a

t



oraz y(a)=a (10)

Na rozwiązanie powyższego zagadnienia będziemy obliczać przybliżone wartości funkcji
yi=y(ti), gdzie ti=a+ih, h=(b-a)/N, dla którego i=(0,1, ..., N), gdzie ti nazywane są punktami
węzłowymi, natomiast h odległością między nimi.
Rozwiniemy y(t) w szereg Taylora w celu wyprowadzenia metody Eulera. Zakładając, że y(t)
jest jedynym rozwiązaniem (10) oraz posiada drugą pochodną, która jest ciągła w przedziale
[a, b] wówczas dla każdego i=1, 2, ..., N.
Wykorzystując powyższe założenia możemy zapisać:

)

(

"

2

)

(

)

(

'

)

(

)

(

)

(

2

1

1

1

i

i

i

i

i

i

i

i

y

t

t

t

y

t

t

t

y

t

y



















(11), gdzie

)

,

(

1





i

i

i

t

t



Oznaczamy h=t

i+1

-t

i

, wówczas otrzymujemy:

)

(

"

2

)

(

'

)

(

)

(

2

1

i

i

i

i

y

h

t

hy

t

y

t

y











(12)

Użyjemy podstawienia y’(t)=f(t, y), wówczas otrzymujemy:

)

(

"

2

))

(

,

(

)

(

)

(

2

1

i

i

i

i

i

y

h

t

y

t

hf

t

y

t

y











(13)

Zapisując ω

i

≈y(t

i

) oraz pomijając błąd przybliżenia otrzymujemy ω

i+1

=ω

i

+hf(t

i

, ω

i

) (14)

76

Powyższy wzór nazywamy metodą Eulera – wzór ten nazywany jest inaczej równaniem
różniczkowym, gdyż można zapisać:

h

t

f

i

i

i

i













1

)

,

(

(15)

Aby wyznaczyć wartość szukanej funkcji y(x) w następnym kroku h, wykorzystujemy
poprzednią wartość funkcji oraz wielkości zmian funkcji – dzięki pochodnej. Natomiast
uwzględniając błąd przybliżenia wzór (13) przyjmuje postać:

i

i

i

i

i

h

t

y

t

f















1

))

(

,

(

, gdzie

)

(

"

2

2

i

i

y

h







(16)

Lokalny błąd dyskretyzacji τ

i+1

(h)

)

(

1

)

,

(

(

1

)

(

1

1

1

1



















i

i

i

i

i

i

i

y

h

t

f

y

y

h

h







(17)

dodatkowo możemy określić krok h, dla którego błąd lokalny jest mniejszy od zadanej
dokładności δ

M

2

h





, gdzie

)

(

"

y

max

M

)

t

,

t

(

1

i

i











Globalny błąd dyskretyzacji g(x)
g(t)= ω(t)-y(t) (19)
Dla wybranego punktu t

i

możemy zapisać:

g

i

=ω

i

-y(t

i

) (20), wówczas





1

2

)

(

)

(

0











t

t

L

i

i

i

i

e

L

hM

t

y

g



, gdzie L- liczba Lipschnitz`a

77

Lokalny błąd dyskretyzacji

Globalny błąd dyskretyzacji

nie będę wskazywał na tego, kto nie narysował go ;)))

Stosując metodę Eulera rozwiązać następujące równanie różniczkowe
y’=y-t

2

+1 0<t<2 y(0)=0,5

dla N=10, wtedy h=0,2; t

i

=ih, ω=0,5 oraz ω

i+1

=ω

i

+hf(t

i

, ω

i

), gdzie f(t,y)=y-t

2

+1

ω

i+1

=ω

i

+h(ω

i

-t

2

i

+1)= ω

i

+0,2[ω

i

-0,04i

2

+1]=1,2ω

i

-0,008i

2

+0,02 dla i=0, 1, ..., 9

Rozwiązanie dokładne ma postać
y(t)=(t

2

-1)-0,5e

t

78

Metoda Eulera – wstecz
Powyższe rozważania dotyczyły metody Eulera w przód, ponieważ krok spełniał warunek
h>0. W analogiczny sposób można wyprowadzić metodą Eulera wstecz przyjmując h<0,
wówczas otrzymujemy:

w

i+1

= w

i

+ hf (t

i+1

, w

i+1

)

(22)

w

i+1

(k)

= w

i

+ hf (t

i+1

, w

i+1

(k-1)

)

(23)

Metoda wstecz różni się w stosunku do metody w przód argumentami funkcji f.

Metoda w przód wykorzystuje do obliczenia przybliżenia wartości z poprzedniego kroku,
natomiast metoda wstecz jest równaniem uwikłanym, ponieważ aby obliczyć kolejne
przybliżenie w

i+1

wykorzystujemy wartości z poprzedniego kroku oraz poszukiwaną wartość

w

i+1

. Takiego równania nie można rozwiązać w sposób bezpośredni. Aby rozwiązać takie

równanie (23) należy zastosować proces iteracyjny, czyli poszukujemy kilkakrotnie w

i+1

(k)

,

stojącej po prawej stronie równania podstawiając jako w

i+1

(k-1)

- lewa strona równania, wynik

przybliżenia z poprzedniej iteracji (k-1). Proces trwa do momentu, kiedy spełniony zostanie
warunek:

|w

i+1

(k)

- w

i+1

(k-1)

| ≤ ε (24)

gdzie ε - tolerancja obliczeń

Kody źródłowe:
Metoda Eulera w przód:

Function Euler_order_One_Forward (t0, tn, y0:extended;
n:integer): extended;
var i:integer;
h, wi, wi_1,t :extended;
Begin
h:= (tn-t0)/n;
yi_1:=y0;
for i:=0 to n-1 do
begin
t:=t0+i*h;
wi:=wi_1+h*f(t,wi_1);
wi_1:=wi;
end;

Result:=wi;

end;

79

Przykładowa funkcja definiująca równanie różniczkowe

)

,

( y

t

f

dt

dy



function f(t,y:extended):extended;
begin
result:=y-sqr(t)+1;
end;


Metoda Eulera “wstecz”

function Euler_Order_One_Backward
t0, tn, y0, eps : extended; n:integer):extended;
var i:integer;
h,wi,wi_p,wi_1,t :extended
begin
h:=(tn-t0)/n
wi_1:=y0;
for i:=0 to n-1 do

begin

t:=t0 + i*h;
wi:=wi_1 + h*f (t, wi_1);
repeat
wi_p:=wi;
wi:=wi_1 + h*f (t, wi_p);
until abs(wi-wi_p)<eps;
wi_1:=wi;
end;
Result:=wi;
end;

Przykładowa funkcja definiująca równanie różniczkowe

)

,

( y

t

f

dt

dy



function f(t,y:extended):extended;
begin
result:=y-sqr(t)+1;
end;


Metody Rungego-Kutty
Powszechnie na całym świecie stosuje się metody Rungego–Kutty czwartego rzędu. Polegają
one na rozwiązaniu zagadnienia:

)

,

( y

t

f

dt

dy



gdzie t





a,b



oraz y(a)=



(25)

i przedstawieniu różnicy funkcji y(t) w punktach t

i+1

oraz t

i

w postaci:

w

i+1

- w

i

=





m

i

j

j

k

c

1

lub inaczej w

i+1

– w

i

= h0 (t

i

,w

i

,h) (26)

80

gdzie m oznacza rząd metody, c

j

są stałymi, a





1

1

1

,

,

,











































j

j

i

j

i

j

i

l

li

i

j

i

j

k

w

t

hf

k

w

h

t

hf

k

gdzie α

j

, β

j

, γ

j

,

δlj

– stałe

h - krok całkowania

Określenie stałych c

j

, α

j

, β

j

otrzymujemy poprzez rozwinięcie funkcji f(t,y) w szereg Taylora

w otoczeniu punktu t

i

Metody R-K – metoda 2 rzędu

Metoda wyprowadzona przez rozwinięcie f(t,y) w szereg Taylora 2 rzędu pozwala określić
stałe c

1

,



1

,



1

, c

2

,



2

,



2

:

Metoda punktu środkowego:

C

1

= 0

C

2

= 1

α

1

= 0

α

2

= h/2



1

= 0



2

= ½ (28)

wówczas możemy zapisać:
k

1

=hf (t

i

,w

i

) k

2

=hf (t

i

+ ½ h, w

i

+ ½ k

1

)

ostatecznie:
w

i+1

=w

i

+k

2

lub
w

i+1

=w

i

+ hf(t

i

+ h/2, w

i

+ h/2 * f(t

i

, w

i

))

81

Interpretacja graficzna:
f

1

=f(t

i

,w

i

) f

2

=f(t

i

+ h/2, w

i

+ h/2 * f

1

)

Stosując metody Rungego–Kutty 2 rzędu , rozwiązać następujące równanie różniczkowe:

y’=y-t

2

+ 1 0



t



2 y(0)=0,5

dla N = 10 wtedy h = 0,2; t

i

= i



h ; w

0

= 0,5 oraz

metoda punktu środkowego w

i+1

= 1,22 w

i

– 0,0088 i

2

+ 0,218

Metoda zmodyfikowana Eulera
w

i+1

= 1,22 w

i

- 0,0088 i

2

+ 0,216

Metoda Heana dla i = 0, 1, ..., 9
w

i+1

= 1,22 w

i

- 0,0088 i

2

+ 0,2173

Rozwiązanie dokładne ma postać:
y(t) = (t

2

+1) - 0,5 e

t

Metody R-K rzędu 4

Metoda wyprowadzona przez rozwinięcie f(t,y) w szereg Taylora 4 rzędu, pozwala określić
stałe we wzorze (wzór ogólny R-K). Poniżej przedstawiono najczęściej stosowaną metodę 4
rzędu

k

1

= hf (t

i

, w

i

)

k

2

= hf (t

i

+ ½ *h, w

i

+ ½ *k

1

)

k

3

= hf (t

i

+ ½ *h, w

i

+ ½ *k

2

) (37)

82

k

4

= hf (t

i

+h, w

i

+k

3

)

ostatecznie:
w

i+1

=w

i

+ 1/6 (k

1

+ 2k

2

+ 2k

3

+ k

4

) (38)

Interpretacja graficzna:

f

1

= f (t

i

, w

i

)

f

2

= f (t

i

+ ½ * h, w

i+1

+ h f

1

)

f

3

= f (t

i

+ ½ * h, w

i

+ ½ * h f

2

)

f

4

= f (t

i

+ h, w

i

+ h f

3

)

linia 3-4 nie jest łamana – jest prosta !!! (na 99%) :)

Metoda R-K jest najczęściej stosowaną metodą do rozwiązywania układów równań
różniczkowych

Przykład:
Stosując metodę R-K 4-tego rzędu rozwiązać następujące równanie różniczkowe:
y’ = y-t

2

+ 1 0



t



2 y(0) = 0,5

dla N = 10 wtedy h = 0,2 t

1

= i*h w

0

= 0,5 oraz:

Metoda Rungego-Kutty 4-rzedu:

w

i+1

= 1,2214 w

i

- 0,008856 i

2

+ 0,00856 i + 0,218593

dla i = 0, 1, ..., 9

Rozwiązanie dokładne ma postać y(t) = (t

2

+1) - 0,5e

t

83

Kody źródłowe:

function RK_Rzedu_4 (t0, tn, y0: extended;
n:integer):extended;
var i:integer;
h, yi, y_i, t, k1, k2, k3, k4: extended;
begin
h:= (tn-t0)/n
y_i = y0;
for i:=0 to n-1 do
begin
t:=t0 + i*h;

k

1 :

= h*ft, y_i);

k

2 :

= h*f(t + 0.5 * h, y_i+0.5 * k1);

k

3 :

= h*f(t+0.5*h, y_i + 0.5 * k

2

);

k

4 :

= h*f(t+h, y_i + k3);

yi:= y_i+(k1 + 2*k2 + 2*k3 + k4)/6;
y_i:=yi;
end;
Result := yi;
end;

Metody wykorzystujące więcej niż jedną przybliżoną wartość dla poprzednich punktów
węzłowych w celu przybliżenia wartości w następnym punkcie węzłowym nazywamy
metodami wielokrokowymi.

Ogólny wzór na m-krotną metodę wielokrokową dla rozwiązania równania

)

,

( y

t

f

dt

dy



gdzie t





a,b



oraz

y(a) =



(55)

jest równaniem różniczkowym dla znalezienia przybliżenia w

i+1

, dla punktu węzłowego t

i+1

,

gdzie m>1, które przyjmuje postać:
w

i+1

= a

m-1

* w

i

+ a

m-2

* w

i-1

+ ... + a

0

* w

i+1_m

+ h [b

n

f(t

i+1

,w

i+1

) + b

n-1

f(t

i

, w

i

) + ... + t + b

0

f(t

i+1_m

, w

i+1_m

)]

dla i = (m-1, m, ..., N-1),
N - liczba węzłów w przedziale, gdzie h = (b-a) / N

84

Równania różniczkowe cząstkowe z warunkami brzegowymi – rodzaje:

a)

równania różniczkowe eliptyczne

•

równanie Laplace’a,

•

równanie Poisson’a,

b)

równanie różniczkowe paraboliczne

•

równanie dyfuzji

c) równanie różniczkowe hiperboliczne
•

równanie falowe

Równania różniczkowe cząstkowe z warunkiem brzegowym – wstęp

Równania różniczkowe prezentowane w poprzednich rozdziałach pozwalały na wyznaczenie
poszukiwanej funkcji jednej zmiennej dla podanych warunków początkowych lub
brzegowych. Problemy przedstawione w tym rozdziale będą opisywane funkcjami wielu
zmiennych.
W zależności od specyfiki problemu można opisać go za pomocą odpowiedniego
równania różniczkowego: Laplace’a, Poisson’a, Dyfuzji, Falowego, itp… . Aby rozwiązać
problem opisany odpowiednim równaniem różniczkowym należy określić warunki brzegowe.
Warunki brzegowe mogą być opisane w formie warunków Dirichleta lub Neumanna.

Równania różniczkowe cząstkowe – warunek brzegowy Dirichleta

Jeżeli funkcja u(x,y) jest poszukiwana na obszarze R, wówczas musi mieć zdefiniowane
warunki w postaci funkcji g(x,y) na brzegu S.
Warunek brzegowy Dirichleta

   

y

,

x

g

y

,

x

u



dla

 

S

y

,

x



gdzie: u(x,y) – poszukiwana funkcja w punktach wewnątrz obszaru R, g(x,y) – zadana
funkcja dla punktów (x,y) należących do brzegu S obszaru R.
Wynika z tego, że w trakcie obliczeń wartości funkcji g w poszczególnych punktach (x,y)



S

nie mogą się zmienić.

85

Równania różniczkowe cząstkowe – warunek brzegowy Neumanna
Warunek brzegowy Neumanna

)

y

,

x

(

g

)

y

,

x

(

n

u







dla

 

S

y

,

x



gdzie:

)

y

,

x

(

n

u





– pochodna normalna poszukiwanej funkcji w punktach należących do

brzegu S obszaru R, g(x,y) – zadana funkcja dla punktów (x,y) należących do brzegu S
obszaru R.

Oznacza to, że na brzegu obszaru zmienność funkcji u(x,y) w kierunku normalnym dla
punktów (x,y)



S jest równa funkcji g(x,y).

Dla zmiennej x

)

y

,

x

(

g

h

)

y

,

h

x

(

u

)

y

,

x

(

u







Dla zmiennej y

)

y

,

x

(

g

h

)

h

y

,

x

(

u

)

y

,

x

(

u







Równania różniczkowe cząstkowe z warunkiem brzegowym - Równania eliptyczne -
Poissona i Laplace’a

Równania różniczkowe cząstkowe eliptyczne znane jako równanie Poissona , dla dwóch
wymiarów i prostokątnego układu współrzędnych przyjmuje postać:

)

y

,

x

(

f

)

y

,

x

(

dy

u

d

)

y

,

x

(

dx

u

d

)

y

,

x

(

u

2

2

2

2

2









(1)

na





d

y

c

b

x

a

y

x

R











;

)

,

(

z warunkami brzegowymi

)

,

(

)

,

(

y

x

g

y

x

u



dla

S

y

x



)

,

(

W powyższym równaniu funkcja f opisuje dane wejściowe na płaszczyźnie dla obszaru R
ograniczonego brzegiem S. Równania tego typu są stosowane dla różnych problemów
fizycznych, które są niezależne od czasu. Najczęściej stosowane są do obliczenia rozkładu
potencjału (np. temperatury) w stanie ustalonym.

Szczególnym przypadkiem równania Poissona gdy f(x,y)=0 jest równanie Laplace’a.

Aby rozwiązać równanie cząstkowe eliptyczne (1) zastosujemy metodę różnic skończonych
MRS.
W tym celu weźmiemy dwie liczby całkowite m i n, i zdefiniujemy kroki całkowania h = (b-
a)/n oraz k = (c-d)/m. Dzięki temu możemy podzielić: przedział [a,b] na n równych części o
szerokości h oraz przedział [c,d] na m równych części o szerokości k.
Umieśćmy siatkę na obszarze R poprzez narysowanie pionowych i poziomych linii
przechodzących przez punkty (xi, yj), takich że:

86

ih

a

x

i





dla i = 0,1, .. n

oraz

jk

c

y

j





dla j = 0,1, .. m

Linie x = xi oraz y = yj nazywamy liniami siatki, natomiast ich punkty przecięcia (xi,yj)
nazywamy punkami węzłowymi siatki. Dla wewnętrznych punktów węzłowych (xi,yj),
które nie należą do brzegu obszaru R, zastosujemy metodę różnic skończonych – MRS.
Metoda ta opiera się na zastąpieniu pochodnych cząstkowych rozwinięciem funkcji w szereg
Taylora w otoczeniu punktu (xi,yj) – (patrz druga pochodna numeryczna) . Wówczas
otrzymujemy:

Dla zmiennej x

)

,

(

12

)

,

(

)

,

(

2

)

,

(

)

,

(

4

4

2

2

1

1

2

2

j

i

j

i

j

i

j

i

j

i

y

x

u

h

h

y

x

u

y

x

u

y

x

u

y

x

x

u























(2)

gdzie:

)

,

(

1

1







i

i

i

x

x



Dla zmiennej y

)

,

(

12

)

,

(

)

,

(

2

)

,

(

)

,

(

4

4

2

2

1

1

2

2

j

i

j

i

j

i

j

i

j

i

x

y

u

k

k

y

x

u

y

x

u

y

x

u

y

x

y

u























(3)

gdzie:

)

,

(

1

1







j

j

j

y

y



87

Podstawiając wzory (2) oraz (3) do równania Poissona (1) otrzymujemy:

2

1

1

2

1

1

)

,

(

)

,

(

2

)

,

(

)

,

(

)

,

(

2

)

,

(

k

y

x

u

y

x

u

y

x

u

h

y

x

u

y

x

u

y

x

u

j

i

j

i

j

i

j

i

j

i

j

i



















)

,

(

12

)

,

(

12

)

,

(

4

4

2

4

4

2

j

i

j

i

x

y

u

k

y

x

u

h

y

x

f



















(4)

dla i = 1,2 … n-1 oraz j = 1,2, … m-1
natomiast warunki brzegowe mają postać:

)

,

(

)

,

(

0

0

j

j

y

x

g

y

x

u



oraz

)

,

(

)

,

(

j

n

j

n

y

x

g

y

x

u



dla j = 0,1, .. m (5)

)

,

(

)

,

(

m

i

m

i

y

x

g

y

x

u



oraz

)

,

(

)

,

(

0

0

y

x

g

y

x

u

i

i



dla i = 1,2, .. n (6)

Ostatecznie wzór na MRS, podstawiając za (xi,yj) = wij oraz wyłączając oszacowanie błędu
zapisujemy:

)

,

(

)

(

)

(

1

2

2

1

,

1

,

2

,

1

,

1

2

j

i

j

i

j

i

j

i

j

i

ij

y

x

f

h

w

w

k

h

w

w

w

k

h































































(7)

dla i = 1,2 … n-1 oraz j = 1,2, … m-1

natomiast warunki brzegowe mają postać:

)

,

(

0

0

j

j

y

x

g

w



oraz

)

,

(

j

n

nj

y

x

g

w



dla j = 0,1, .. m (8)

)

,

(

0

0

y

x

g

w

i

i



oraz

)

,

(

m

i

im

y

x

g

w



dla i = 1,2, .. n (9)

gdzie w

ij

jest przybliżeniem u(x

i

,y

j

).

88

Analizując równanie (7) można zauważyć, że w celu wyznaczenia przybliżenie rozwiązania w
punkcie (xi,yj), potrzebne są wartości przybliżenia rozwiązania w czterech sąsiednich
punktach – patrz rysunek obok.

Postać macierzowa
Jak widać poszukujemy rozwiązania równania Poissona dla wewnętrznych węzłów siatki, w
tym celu dla każdego punktu siatki układamy równanie (7). Po uwzględnieniu warunków
brzegowych otrzymujemy układ (n-1)



(m-1) równań liniowych. Układ równań możemy

rozwiązać metodami deterministycznymi bądź iteracyjnymi.

Zakładając, że h = k możemy wyprowadzić uproszczoną postać równania (7) :

)

,

(

)

(

)

(

4

2

1

,

1

,

,

1

,

1

j

i

j

i

j

i

j

i

j

i

ij

y

x

f

h

w

w

w

w

w





















(10)

dla i = 1,2 … n-1 oraz j = 1,2, … m-1


natomiast warunki brzegowe są takie same jak (8) i (9) :

Postać iteracyjna

Równanie (7) i (10) można przedstawić również w postaci iteracyjnej, wyprowadzając wij z
każdego równania, gdzie l oznacza numer iteracji :

z równania (7)

 

 





1

2

)

(

)

(

)

,

(

2

)

(

1

,

)

(

1

,

2

)

(

,

1

)

(

,

1

2

)

1

(

























k

h

l

j

i

l

j

i

k

h

l

j

i

l

j

i

j

i

l

ij

w

w

w

w

y

x

f

h

w

(11)

z równania (10)

4

)

(

)

,

(

)

(

1

,

)

(

1

,

)

(

,

1

)

(

,

1

2

)

1

(

l

j

i

l

j

i

l

j

i

l

j

i

j

i

l

ij

w

w

w

w

y

x

f

h

w























(12)

89

Równania eliptyczne - przykład 1

Wyznaczyć rozkład temperatury w stanie ustalonym dla cienkiej kwadratowej
metalowej płytki o wymiarach 0,5 m na 0,5 m. Na brzegu płytki znajdują się źródła
ciepła utrzymujące temperaturę na poziomie 0°C dla boku dolnego i prawego, natomiast
temperatura boku górnego i lewego zmienia się liniowo od 0°C do 100°C Problem
rozwiązać układając układ równań liniowych (postać macierzowa) dla wewnętrznych
węzłów siatki 5 x 5 – układ równań rozwiązać metodą Gaussa -Siedla.

Siatka dyskretyzacyjna 5 x 5

Problem ten opisuje równanie Laplace’a

0

)

,

(

)

,

(

)

,

(

2

2

2

2

2









y

x

dy

T

d

y

x

dx

T

d

y

x

T

z warunkami brzegowymi:
1) T(0,y) = 0 [°C ]
2) T(x,0) = 0 [°C ]
3) T(0,5,y) = 200y [°C ]
4) T(x,0,5) = 200x [°C ]

Układ równań tworzymy w oparciu o wzór (10)



























































































































































































































25

0

0

50

0

0

100

50

25

4

1

0

1

0

0

0

0

0

1

4

1

0

1

0

0

0

0

0

1

4

0

0

1

0

0

0

1

0

0

4

1

0

1

0

0

0

1

0

1

4

1

0

1

0

0

0

1

0

1

4

0

0

1

0

0

0

1

0

0

4

1

0

0

0

0

0

1

0

1

4

1

0

0

0

0

0

1

0

1

4

3

,

1

2

,

1

1

,

1

3

,

2

2

,

2

1

,

2

3

,

3

2

,

3

1

,

3

T

T

T

T

T

T

T

T

T

90

Rozwiązanie układu równań metodą Gaussa-Siedla, daje wyniki

75

,

18

50

,

12

25

,

6

50

,

37

00

,

25

50

,

12

25

,

56

50

,

38

75

,

18

3

,

1

2

,

1

1

,

1

3

,

2

2

,

2

1

,

2

3

,

3

2

,

3

1

,

3

T

T

T

T

T

T

T

T

T

Ostatecznie wyniki dla założonej siatki mają postać:

00

,

0

00

,

0

00

,

0

00

,

0

00

,

0

0

00

,

25

75

,

18

50

,

12

25

,

6

00

,

0

125

,

0

00

,

50

50

,

37

00

,

25

50

,

12

00

,

0

25

,

0

00

,

75

25

,

56

50

,

37

75

,

18

00

,

0

375

,

0

00

,

100

00

,

75

00

,

50

00

,

25

00

,

0

5

,

0

5

,

0

375

,

0

25

,

0

125

,

0

0





















y

y

y

y

y

x

x

x

x

x

Rozwiązanie dokładne problemu ma postać:

xy

y

x

T

400

)

,

(



Natomiast oszacowanie błędu obliczeń :

0

)

,

(

)

,

(

12

4

4

4

4

2



























j

i

j

i

x

y

T

y

x

T

h







, bo

0

)

400

(

)

,

(

4

4

4

4













xy

x

y

x

x

T

i

0

)

400

(

)

,

(

4

4

4

4













xy

y

y

x

y

T

Wynika z tego, że rozwiązanie numeryczne jest równe z rozwiązaniu dokładnemu.

Równania eliptyczne - przykład 1 - kody źródłowe

Funkcje charakteryzujące dane równanie różniczkowe eliptyczne

function F( X,Y : real ) : real;
begin
F := 0;
end;

function G( X,Y : real ) : real;
begin
if X = 0 then G := 0;
if X = 0,5 then G := 200*Y;
if Y = 0 then G := 0;
if Y = 0,5 then G := 200*X;
end;

Definicja wymaganych stałych i typów

91

const max = 4;
type
Macierz = array [ 0..max, 0..max ] of real;
Vector = array [ 0..max ] of real;


Funkcja obliczająca przybliżenie równania różniczkowego eliptycznego

function PoissonEq(A,B,C,D,TOL : real;N,M,NMAX : integer; var W : Macierz)
:boolean;
var
X,Y : Vector; H,K,V,VV,Z,E : real;
M1,M2,N1,N2,I,J,L,LL: integer;
Ok :boolean;
begin
M1 := M - 1; M2 := M - 2;
N1 := N - 1; N2 := N - 2;
H := ( B - A ) / N; K := ( D - C ) / M;
for I := 0 to N do X[I] := A + I * H;
for J := 0 to M DO Y[J] := C + J * K;
for I := 1 to N1 do begin
W[I,0] := g(X[I],Y[0]);
W[I,M] := g(X[I],Y[M])
end;
for J := 0 to M do begin
W[0,J] := g(X[0],Y[J]);
W[N,J] := g(X[N],Y[J])
end;
for I := 1 to N1 do for J := 1 to M1 do W[I,J] := 0.0;
V := H * H / ( K * K ); VV := 2.0 * ( 1.0 + V );
L := 1; OK := false;
while ( ( L <= NMAX ) and ( not OK ) ) do begin
Z := (-H*H*F(X[1],Y[M1])+G(A,Y[M1])+V*G(X[1],D)+V*W[1,M2]+W[2,M1])/VV;
E := abs( W[1,M1] - Z ); W[1,M1] := Z;
for I := 2 to N2 do begin
Z := (-H*H*F(X[I],Y[M1])+V*G(X[I],D)+W[I-1,M1]+W[I+1,M1]+V*W[I,M2])/VV;
if ( abs( W[I,M1] - Z ) > E ) then E := abs( W[I,M1] - Z );
W[I,M1] := Z
end;
Z := (-
H*H*F(X[N1],Y[M1])+G(B,Y[M1])+V*G(X[N1],D)+W[N2,M1]+V*W[N1,M2])/VV;
if ( abs( W[N1,M1] - Z ) > E ) then E := abs( W[N1,M1] - Z ); W[N1,M1] := Z;
for LL := 2 to M2 do begin
J := M2 - LL + 2; Z := (-H*H*F(X[1],Y[J])+G(A,Y[J])+V*W[1,J+1]+V*W[1,J-
1]+W[2,J])/VV;
if ( abs( W[1,J] - Z ) > E ) then E := abs( W[1,J] - Z );
W[1,J] := Z;
for I := 2 to N2 do begin
Z := (-H*H*F(X[I],Y[J])+W[I-1,J]+V*W[I,J+1]+V*W[I,J-1]+W[I+1,J])/VV;
if ( abs( W[I,J] - Z ) > E ) then E := abs( W[I,J] - Z );

92

W[I,J] := Z
end;
Z := (-H*H*F(X[N1],Y[J])+G(B,Y[J])+W[N2,J]+V*W[N1,J+1]+V*W[N1,J-1])/VV;
if ( abs( W[N1,J] - Z ) > E ) then E := abs( W[N1,J] - Z ); W[N1,J] := Z
end;
Z := ( -H * H * F( X[1],Y[1] ) + V * G( X[1], C ) + G( A, Y[1] ) + V * W[1,2] + W[2,1] ) /
VV;
if ( abs( W[1,1] - Z ) > E ) then E := abs( W[1,1] - Z );
W[1,1] := Z;
for I := 2 to N2 do begin
Z := (-H*H*F(X[I],Y[1])+V*G(X[I],C)+ W[I+1,1]+W[I-1,1]+V*W[I,2])/VV;
if ( abs( W[I,1] - Z ) > E ) then E := abs( W[I,1] - Z );
W[I,1] := Z
end;
Z := (-H*H*F(X[N1],Y[1])+V*G(X[N1],C)+G(B,Y[1])+W[N2,1]+V*W[N1,2])/VV;
if ( abs( W[N1,1] - Z ) > E ) then E := ABS( W[N1,1] - Z );
W[N1,1] := Z;
if ( E <= TOL ) then OK := true
else L := L + 1;
end;
Result := OK;
end;

93

Równania eliptyczne - przykład 2
Wyznaczyć rozkład temperatury w stanie ustalonym dla cienkiej kwadratowej metalowej
płytki o wymiarach 1 m na 1 m. Na trzech brzegach płytki znajdują się źródła ciepła
utrzymujące temperaturę na poziomie 100°C dla boku górnego, 50°C dla boku lewego, 0°C
dla boku dolnego (warunki Dirichleta), natomiast przy czwartym lewym boku jest
umieszczony izolator termiczny (warunek Neumanna – pochodna normalna równa jest 0).
Problem rozwiązać wg postaci iteracyjnej dla wewnętrznych węzłów siatki 101 x 101.
Tolerancję obliczeń przyjąć na poziomie 10

-6

.

Problem ten opisuje równanie Laplace’a

0

)

,

(

)

,

(

)

,

(

2

2

2

2

2









y

x

dy

T

d

y

x

dx

T

d

y

x

T

z warunkami brzegowymi:
1) T(0,y) = 0 [°C ]
2) T(x,0) = 50 [°C ]
3) T(1,y) = T(1-h,y) [°C ]
4) T(x,1) = 100 [°C ]

94

Przybliżenie rozwiązania równania w postaci mapy barwnej po 13484 iteracjach

Tabela wyników dla wybranych punków prawego boku płytki - izolator

Numer węzła

Numer iteracji

Dokładne

0

1000

3000

7000

13000

T100,0

100,00

100,00

100,00

100,00

100,00

100,00

T100,10

0,00

75,17

86,03

89,59

89,92

90,00

T100,20

0,00

52,58

72,46

79,24

79,87

80,00

T100,30

0,00

34,03

59,66

69,00

69,86

70,00

T100,40

0,00

20,27

47,91

58,89

59,90

60,00

T100,50

0,00

11,08

37,37

48,90

49,97

50,00

T100,60

0,00

5,55

28,07

39,03

40,05

40,00

T100,70

0,00

2,55

19,91

29,22

30,09

30,00

T100,80

0,00

1,07

12,70

19,47

20,10

20,00

T100,90

0,00

0,37

6,17

9,73

10,06

10,00

T100,100

0,00

0,00

0,00

0,00

0,00

0,00

95

Równania eliptyczne - przykład 2 - kody źródłowe


Funkcje charakteryzujące dane równanie różniczkowe eliptyczne

function F( X,Y : real ) : real;
begin
F := 0;
end;


Definicja wymaganych stałych i typów

const max = 100;
type
Macierz = array [ 0..max, 0..max ] of real;
Vector = array [ 0..max ] of real;

function PoissonEquationIteration (a,b,c,d, TOL : real ; N,M ,nmax : integer; var V :
Macierz) : boolean;
var EPS,H,K,Vs : real; X,Y : Vector; I,J,L : integer;
begin
H := ( B - A ) / N; K := ( D - C ) / M;
for I := 0 to N do X[I] := A + I * H; for J := 0 to M do Y[J] := C + J * K;
for i:=0 to N do for j:=0 to M do V[i,j]:=0;
L := 0; Result := true;
repeat
eps:=0; L := L + 1;
for i:=0 to N do for j:=0 to M do begin
Vs:=V[i,j];
if j=0 then V[i,j]:=0 else // bok dolny
if j=M then V[i,j]:=100 else // bok gora
if i=0 then V[i,j]:= 50 else // bok lewo
if i=N then V[i,j]:= V[i-1,j] else // bok prawo
V[i,j]:=(- h*h* F(X[i],Y[j] ) + (V[i+1,j]+V[i-1,j])+sqr(h/k)*(V[i,j+1]+V[i,j-1])
)/(2*(sqr(h/k)+1));
If abs(V[i,j]-Vs) > eps then eps := abs(V[i,j]-Vs);
end;
until (eps< TOL ) or ( L > nmax);
if (L > nmax) then Result := false;
end;

96

Równania różniczkowe cząstkowe z warunkiem brzegowym

Równania paraboliczne - Dyfuzja

Równania różniczkowe cząstkowe paraboliczne znane jako równanie dyfuzji lub
przewodnictwa, w zależności od jednego wymiaru oraz czasu przyjmuje postać:

dla oraz

(13)

z warunkami brzegowymi dla


i początkowymi dla


Równania tego typu są stosowane dla różnych problemów fizycznych, które są zależne od
czasu. Najczęściej stosowane są do obliczenia przepływu ciepła wzdłuż pręta o długości l,
przy założeniu że, powierzchnia boczna pręta jest odizolowana od otoczenia. Stała



jest

niezależna od położenia oraz czasu i najczęściej określa przewodność cieplną materiału z
którego zrobiony jest pręt. Funkcja f określa początkowy rozkład temperatury w pręcie.

Aby rozwiązać równanie cząstkowe paraboliczne (13) zastosujemy metodę różnic
skończonych MRS.

W tym celu weźmiemy liczbę
całkowitą m i zdefiniujemy
krok całkowania h = (b-a)/m.
Dodatkowo ustalimy wartość
kroku czasowego k. Dzięki
temu możemy wyznaczyć
węzły siatki (xi , tj ), gdzie:

dla i = 0,1, .. m oraz dla j = 0,1, ..


Dla wewnętrznych punktów węzłowych (xi,tj), które nie należą do brzegu, zastosujemy
metodę różnic skończonych – MRS. Metoda ta opiera się na zastąpieniu pochodnych
cząstkowych rozwinięciem funkcji w szereg Taylora w otoczeniu punktu (xi,tj) – (patrz
pierwsza i druga pochodna numeryczna) . Wówczas otrzymujemy:

Dla zmiennej x (14)

gdzie:

Dla zmiennej t (15)

gdzie:

0

)

,

(

)

,

(

2

2

2





t

x

dx

u

d

t

x

dt

du



l

x





0

0



t

0

)

,

(

)

,

0

(





t

l

u

t

u

0



t

)

(

)

0

,

(

x

f

x

u



l

x





0

ih

x

i



jk

t

j



)

,

(

12

)

,

(

)

,

(

2

)

,

(

)

,

(

4

4

2

2

1

1

2

2

j

i

j

i

j

i

j

i

j

i

t

x

u

h

h

t

x

u

t

x

u

t

x

u

y

x

x

u























)

,

(

1

1







i

i

i

x

x



)

,

(

2

)

,

(

)

,

(

)

,

(

2

2

1

j

i

j

i

j

i

j

i

x

t

u

k

k

t

x

u

t

x

u

t

x

t

u



















)

,

(

1





j

j

j

t

t



97





j

i

j

i

j

i

j

i

w

w

h

k

w

h

k

w

,

1

,

1

2

2

,

2

2

1

,

2

1































Podstawiając wzory (14) oraz (15) do równania dyfuzji (13) oraz wyłączając oszacowanie
błędu zapisujemy:

(16)

warunek brzegowy : (17)

dla i = 1,2 … m-1 oraz j = 1,2, …

natomiast lokalny błąd odcięcia jest równy :

(18)


gdzie: oraz

Ostatecznie wzór na MRS, podstawiając za ( xi , tj ) = wi,j oraz wyłączając oszacowanie
błędu zapisujemy:

(19)


przekształcając powyższy wzór ze względu na wi,j+1 ,wówczas otrzymujemy:


(20)



lub

gdzie: (21)

dla i = 1,2 … m-1 oraz j = 1,2, …

0

)

,

(

)

,

(

2

)

,

(

)

,

(

)

,

(

2

1

1

2

1

















h

t

x

u

t

x

u

t

x

u

k

t

x

u

t

x

u

j

i

j

i

j

i

j

i

j

i



0

)

,

(

)

,

0

(





j

j

t

l

u

t

u

)

,

(

12

)

,

(

2

4

4

2

2

2

2

j

i

j

i

ij

t

x

u

h

x

t

u

k





















)

,

(

1

1







i

i

i

x

x



)

,

(

1





j

j

j

t

t



0

2

2

,

1

,

,

1

2

,

1

,

















h

w

w

w

k

w

w

j

i

j

i

j

i

j

i

j

i











j

i

j

i

j

i

j

i

w

w

w

w

,

1

,

1

,

1

,

2

1



















2















h

k





98

Równania paraboliczne – Schemat jawny

Analizując równanie (21) można zauważyć,
że w celu wyznaczenia przybliżenia
rozwiązania w punkcie (xi,tj+1),
potrzebne są wartości trzech punktów z
poprzedniego kroku czasowego tj - patrz
rysunek obok. Metoda ta nazywana jest
rozwiązywaniem równania parabolicznego
wg schematu jawnego


Jak

widać

poszukujemy

rozwiązania

równania dyfuzji dla wewnętrznych węzłów siatki, w tym celu dla każdego punktu siatki
układamy równanie (21). Obliczenia rozpoczynamy od kroku czasowego j = 1, gdyż najpierw
musimy wygenerować wartości początkowe dla t = 0 dla wszystkich punktów u( xi , 0).
Obliczmy je zgodnie z warunkami początkowymi u( xi , 0) = f(xi).

A więc mamy do rozwiązania układ (m-1) równań, który powstał poprzez ułożenie równania
(21) dla każdego punktu wi,j siatki dla i = 1,2, .. m-1 :

(22)

co rozpisując daje nam:



(23)





Czyli w ten sposób możemy wyznaczyć poszukiwane wartości w(j+1). Wystarczy tylko
przemnożyć macierz A przez wektor w(j) - wartości dla poprzedniego kroku.

)

(

)

1

(

j

j

Aw

w





























































































































j

m

j

j

j

m

j

j

w

w

w

w

w

w

,

1

,

2

,

1

1

,

1

1

,

2

1

,

1

)

2

1

(

0

0

0

0

0

0

)

2

1

(





















99

Równania paraboliczne - Dyfuzja - przykład 1

Rozwiązać poniższe równanie przewodnictwa:

dla oraz


z warunkami brzegowymi dla

i początkowymi dla

Obliczenia wykonać dla czasu t = 0.5 używając schematu jawnego najpierw dla h = 0.1 oraz k
= 0.0005 , następnie dla h = 0.1 oraz k = 0.01 ,przyjąć



2 = 1

Rozwiązanie dokładne ma postać:



Tabela wyników – schemat jawny

Mały krok czasu

Duży krok czasu

xi

u( xi , 0.5 )

wi,1000
k = 0.0005

| u( xi , 0.5 ) –
wi,1000|

wi,50
k = 0.01

| u( xi , 0.5 ) –
wi,50|

0,0

0

0

0

0,1

0,00222241 0,00228652 6,411



10-5

8,19876



107 8,199



107

0,2

0,00422728 0,00434922 1,219



10-4

-1,55719



108 1,557



108

0,3

0,00581836 0,00598619 1,678



10-4

2,13833



108 2,138



108

0,4

0,00683989 0,00703719 1,973



10-4

-2,50642



108 2,506



108

0,5

0,00719788 0,00739934 2,075



10-4

2,62685



107 2,627



108

0,6

0,00683989 0,00703719 1,973



10-4

-2,49015



108 2,490



108

0,7

0,00581836 0,00598619 1,678



10-4

2,11200



107 2,112



108

0,8

0,00422728 0,00434922 1,219



10-4

-1,53086



108 1,531



108

0,9

0,00222241 0,00228652 6,411



10-5

8,03604



107 8,036



107

1,0

0

0

0

0

)

,

(

)

,

(

2

2

2





t

x

dx

u

d

t

x

dt

du



1

0





x

0



t

0

)

,

1

(

)

,

0

(





t

u

t

u

0



t

)

sin(

)

0

,

(

x

x

u





1

0





x

)

sin(

)

,

(

2

x

e

t

x

u

t









100

Interpretacja graficzna




























Równania paraboliczne – Schemat jawny

Analizując wyniki z przykładu 1 można zauważyć, że metoda wykorzystująca schemat jawny
jest warunkowo stabilna, gdyż dla zbyt dużych wartości kroku wyniki są rozbieżne w
stosunku do rozwiązania dokładnego. Wynika to wartości własnej macierzy A , dla której
powinien być spełniony warunek

(23)

z powyższego wyrażenia wynika, że warunek stabilności tej metody jest:

(24)

Wykorzystując powyższy wzór można wyznaczyć dla przykładu 1 graniczną wartość kroku
czasowego, dla której schemat jawny jest stabilny:

wtedy

 





1

sin

4

1

max

)

(

2

2

1

1













m

i

m

i

A







2

1

2

2



h

k



2

2

1

















h

k

0,005

1

1

,

0

2

1

2

















k

101

0

)

,

(

)

,

(

2

)

,

(

)

,

(

)

,

(

2

1

1

2

1

















h

t

x

u

t

x

u

t

x

u

k

t

x

u

t

x

u

j

i

j

i

j

i

j

i

j

i



Równania paraboliczne – Schemat niejawny

Poprzednia metoda wykorzystująca schemat jawny, była warunkowo stabilna, gdyż należy
zwracać uwagę na wielkości kroków całkowania. Metoda opisana poniżej nie ma tej wady, i
jest bezwarunkowo stabilna dla różnych wielkości kroków całkowania. Metoda ta
wykorzystuje schemat niejawny gdyż zalicz się do metod pośrednich wymagających
wielokrotnych iteracji dla danej chwili czasu. Różnica polega na innym określeniu pochodnej
cząstkowej po czasie:

Dla zmiennej t (25)

gdzie:

Wówczas podstawiając wzór (25) oraz (14) do równania parabolicznego (13), pomijając błąd
obcięcia otrzymamy :

(26)

Modyfikując powyższy wzór, podstawiając za ( xi , tj ) = wi,j zapisujemy:

(27)

lub:

gdzie: (28)

dla i = 1,2 … m-1 oraz j = 1,2, …



Analizując równanie (28) można
zauważyć, że w celu wyznaczenia
przybliżenia rozwiązania w punkcie
(xi,tj), potrzebne są wartości trzech
punktów:

dwóch

sąsiednich

z

aktualnego kroku czasowego oraz
wartość dla tego samego węzła, ale
dla poprzedniej iteracji tj-1 - patrz
rysunek obok. Metoda ta nazywana
jest

rozwiązywaniem

równania

parabolicznego

wg

schematu

niejawnego


Jak widać poszukujemy rozwiązania równania dyfuzji dla wewnętrznych węzłów siatki, w
tym celu dla każdego punktu siatki układamy równanie (28). Obliczenia rozpoczynamy od
kroku czasowego j = 1, gdyż najpierw musimy wygenerować wartości początkowe dla t = 0
dla wszystkich punktów u( xi , 0). Obliczmy je zgodnie z warunkami początkowymi u( xi , 0)
= f(xi).

)

,

(

2

)

,

(

)

,

(

)

,

(

2

2

1

j

i

j

i

j

i

j

i

x

t

u

k

k

t

x

u

t

x

u

t

x

t

u



















)

,

(

1

j

j

j

t

t







0

2

2

,

1

,

,

1

2

1

,

,

















h

w

w

w

k

w

w

j

i

j

i

j

i

j

i

j

i



1

,

,

1

,

1

,

)

(

)

2

1

(















j

i

j

i

j

i

j

i

w

w

w

w





2















h

k





102

































































































































1

,

1

1

,

2

1

,

1

,

1

,

2

,

1

)

2

1

(

0

0

0

0

0

0

)

2

1

(

j

m

j

j

j

m

j

j

w

w

w

w

w

w





















A więc mamy do rozwiązania układ (m-1) równań, który powstał poprzez ułożenie równania
(28) dla każdego punktu wi,j siatki dla i = 1,2, .. m-1 :


(29)


co rozpisując daje nam:





(30)







Czyli w ten sposób możemy wyznaczyć poszukiwane wartości w(j). Najlepiej do rozwiązania
tego typu układów użyć rozkładu Crouta dla macierzy trójdiagonalnych


Równania paraboliczne - Dyfuzja - przykład 2 (1)

Rozwiązać poniższe równanie przewodnictwa:

dla oraz


z warunkami brzegowymi dla

i początkowymi dla

Obliczenia wykonać dla czasu t = 0.5 używając schematu niejawnego dla h = 0.1 oraz k
= 0.01 ,przyjąć



2 = 1

Rozwiązanie dokładne ma postać:

)

1

(

)

(





j

j

w

Aw

0

)

,

(

)

,

(

2

2

2





t

x

dx

u

d

t

x

dt

du



1

0





x

0



t

0

)

,

1

(

)

,

0

(





t

u

t

u

0



t

)

sin(

)

0

,

(

x

x

u





1

0





x

)

sin(

)

,

(

2

x

e

t

x

u

t









103

Tabela wyników – schemat niejawny

xi

u( xi , 0.5 )

wi,50
k = 0.01

| u( xi , 0.5 )
– wi,50|

0,0

0

0

0,1

0,00222241 0,00289802 6,756



10-4

0,2

0,00422728 0,00551236 1,285



10-3

0,3

0,00581836 0,00758711 1,769



10-3

0,4

0,00683989 0,00891918 2,079



10-3

0,5

0,00719788 0,00937818 2,186



10-3

0,6

0,00683989 0,00891918 2,079



10-3

0,7

0,00581836 0,00758711 1,769



10-3

0,8

0,00422728 0,00551236 1,285



10-3

0,9

0,00222241 0,00289802 6,756



10-4

1,0

0

0

Równania paraboliczne – Schemat niejawny

Analizując wyniki z przykładu 2 można zauważyć, że metoda wykorzystująca schemat
niejawny jest bezwarunkowo stabilna, gdyż jest ona zawsze zbieżna nie zależnie od wielkości
kroków całkowania. Wynika to ze spełniania poniższego warunku dla promienia spektralnego
macierzy A:

(31)


Równania paraboliczne – kody źródłowe

Funkcja określająca warunki początkowe równania paraboliczne

function F( X,Y : real ) : real;
begin
F := sin( pi * x );
end;

Definicja wymaganych stałych i typów

const
max = 100;
type
Vector = array [ 0..max ] of real;

 





1

sin

4

1

max

)

(

2

2

1

1













m

i

m

i

A







104

Funkcji „HeatingForward” - schemat jawny

function HeatingForward (FX,ALPHA ,FT : real; M,N : integer; var W : Vector):
boolean ;
var
H,K,VV : real; I,J : integer; W_1 : Vector;
begin
H := FX / M; K := FT / N;
VV := ALPHA * ALPHA * K / ( H * H );
if K <= (0.5*SQR(H/ALPHA)) then begin
W_1[0] := 0; W_1[M] := 0; Result := True;
W[0] := 0; W[M] := 0;
for I := 1 to M-1 do W_1[I] := F( I * H );
for J := 1 to N do begin
for I := 1 to M-1 do W[I] := (1.0 - 2.0 * VV)*W_1[I] + VV*( W_1[I-1] + W_1[I+1]
);
W_1 := W;
end;
end
else Result := False;
end;

Funkcji „HeatingBackward” - schemat niejawny

function HeatingBackward (FX,alpha ,FT : real; m,n : integer; var W : Vector):
boolean ;
var
L,U,Z : Vector; H,K,VV : real; I1,I,J : integer;
begin
Result := true;
H := FX / M; K := FT / N;
VV := ALPHA * ALPHA * K / ( H * H );
W[0] := 0; W[M] := 0;

for I := 1 to M-1 do W[I] := F( I * H );

L[1] := 1.0 + 2.0 * VV; U[1] := -VV / L[1];

for I := 2 to M-2 do begin
L[I] := 1.0 + 2.0 * VV + VV * U[I-1]; U[I] := -VV / L[I]
end;

L[M-1] := 1.0 + 2.0 * VV + VV * U[M-2];

for J := 1 to N do begin
Z[1] := W[1] / L[1];

for I := 2 to M-1 do Z[I] := ( W[I] + VV * Z[I-1] ) / L[I];

W[M-1] := Z[M-1];

105

for I1 := 1 to M-2 do begin
I := M-2 - I1 + 1; W[I] := Z[I] - U[I] * W[I+1]
end;

end;
end;

106

Równanie hiperboliczne – Falowe

Równania różniczkowe cząstkowe hiperboliczne znane jako równanie falowe, w zależności
od jednego wymiaru oraz czasu przyjmuje postać:

(32)

dla oraz

z warunkami brzegowymi

dla

i początkowymi

oraz

dla

Równania tego typu są stosowane dla różnych problemów fizycznych, które są zależne od
czasu. Najczęściej stosowane są do obliczenia wartości fali płaskiej w dowolnym punkcie
pomiędzy 0 a l dla dowolnej chwili czasu większej od zera. Funkcje f(x) i g(x) określają
warunki początkowe położenia oraz prędkości rozchodzenia się fali.

Aby rozwiązać równanie cząstkowe paraboliczne (32) zastosujemy metodę różnic
skończonych MRS.

W tym celu weźmiemy liczbę całkowitą m i zdefiniujemy krok całkowania h = (b-a)/m.
Dodatkowo ustalimy wartość kroku czasowego k. Dzięki temu możemy wyznaczyć węzły
siatki (x

i

, t

j

), gdzie:

dla i = 0,1, .. m oraz dla j = 0,1, ..

Dla wewnętrznych punktów węzłowych (x

i

,t

j

), które nie należą do brzegu, zastosujemy

metodę różnic skończonych – MRS. Metoda ta opiera się na zastąpieniu pochodnych
cząstkowych rozwinięciem funkcji w szereg Taylora w otoczeniu punktu (x

i

,t

j

) – (patrz druga

pochodna numeryczna) . Wówczas otrzymujemy:

Dla zmiennej x

(33)

0

)

,

(

)

,

(

2

2

2

2

2





t

x

dx

u

d

t

x

dt

u

d



l

x





0

0



t

0

)

,

(

)

,

0

(





t

l

u

t

u

0



t

)

(

)

0

,

(

x

f

x

u



)

(

)

0

,

(

x

g

x

t

u







l

x





0

ih

x

i



jk

t

j



)

,

(

12

)

,

(

)

,

(

2

)

,

(

)

,

(

4

4

2

2

1

1

2

2

j

i

j

i

j

i

j

i

j

i

t

x

u

h

h

t

x

u

t

x

u

t

x

u

t

x

x

u























107

gdzie

Dla zmiennej t

(34)

gdzie

Podstawiając wzory (33) oraz (34) do równania falowego (31) oraz wyłączając oszacowanie
błędu zapisujemy:

(35)

warunek brzegowy :

(36)

dla

i = 1,2 … m-1 oraz j = 1,2, …

natomiast lokalny błąd odcięcia jest równy :

(37)

gdzie:

oraz

Ostatecznie wzór na MRS, podstawiając za ( x

i

, t

j

) = w

i,j

zapisujemy:

(38)

przekształcając powyższy wzór ze względu na w

i,j+1

,wówczas otrzymujemy:

gdzie:

(39)

dla i = 1,2 … m-1 oraz j = 1,2, …

z warunkami brzegowymi

dla j = 1,2, …

i początkowymi

dla i = 1,2, … m-1

)

,

(

1

1







i

i

i

x

x



)

,

(

12

)

,

(

)

,

(

2

)

,

(

)

,

(

4

4

2

2

1

1

2

2

j

i

j

i

j

i

j

i

j

i

x

t

u

k

k

t

x

u

t

x

u

t

x

u

t

x

y

u























)

,

(

1

1







j

j

j

t

t



0

)

,

(

)

,

(

2

)

,

(

)

,

(

)

,

(

2

)

,

(

2

1

1

2

2

1

1





















h

t

x

u

t

x

u

t

x

u

k

t

x

u

t

x

u

t

x

u

j

i

j

i

j

i

j

i

j

i

j

i



0

)

,

(

)

,

0

(





j

j

t

l

u

t

u

























)

,

(

)

,

(

12

1

4

4

2

2

4

4

2

j

i

j

i

ij

t

x

u

h

x

t

u

k









)

,

(

1

1







i

i

i

x

x



)

,

(

1

1







j

j

j

t

t



0

2

2

2

,

1

,

,

1

2

2

1

,

,

1

,





















h

w

w

w

k

w

w

w

j

i

j

i

j

i

j

i

j

i

j

i











1

,

,

1

,

1

2

,

2

1

,

1

2



















j

i

j

i

j

i

j

i

j

i

w

w

w

w

w





h

k



 

0

,

,

0





j

m

j

w

w

)

(

0

,

i

i

x

f

w



108

Analizując równanie (39) można zauważyć, że w celu wyznaczenia przybliżenia rozwiązania
w punkcie (x

i

,t

j+1

), potrzebne są wartości trzech punktów z poprzedniego kroku

czasowego t

j

oraz jednej wartości dla czasu t

j-1

– patrz rysunek poniżej.

Jak widać poszukujemy rozwiązania równania falowego dla wewnętrznych węzłów siatki, w
tym celu dla każdego punktu siatki układamy równanie (39). Obliczenia rozpoczynamy od
kroku czasowego j = 2, gdyż najpierw musimy wygenerować wartości początkowe dla t = 0
dla wszystkich punktów u( x

i

, 0)=f(x

i

), następnie wygenerujemy z drugiego warunku

początkowego wartości punktów dla t = t

1

jako du( x

i

, t

1

) / dt = g(x

i

).

A więc mamy do rozwiązania układ (m-1) równań, który powstał poprzez ułożenie równania
(39) dla każdego punktu w

i,j

siatki dla i = 1,2, .. m-1 :

(40)

co rozpisując daje nam:

(41)

Czyli w ten sposób możemy wyznaczyć poszukiwane wartości w

(j+1)

. Wystarczy tylko

przemnożyć macierz A przez wektor w

(j)

- wartości dla poprzedniego kroku, a następnie odjąć

od tego wektor wynikowy dla t

j-1

czyli w

(j-1)

.

)

1

(

)

(

)

1

(









j

j

j

w

Aw

w



































































































































































1

,

1

1

,

2

1

,

1

,

1

,

2

,

1

2

2

2

2

2

2

1

,

1

1

,

2

1

,

1

)

1

(

2

0

0

0

0

0

0

)

1

(

2

j

m

j

j

j

m

j

j

j

m

j

j

w

w

w

w

w

w

w

w

w

























109

Ponieważ obliczania rozpoczynamy od kroku czasowego j = 2 musimy oprócz warunków
początkowych dla j = 0 określić wartości dla j = 1, skorzystamy z drugiego warunku
początkowej prędkości fali:

dla

(42)

Zamienimy pochodną cząstkową po czasie wykorzystując wzór schematu jawnego dla
równań parabolicznych:

(43)

Wyprowadźmy u(x

i

,t

1

) z powyższego równania, wtedy otrzymujemy:

(44)

Podstawmy do równania (44) wartość pochodnej z równania (42) :

(45)

Podstawmy za u( x

i

,t

j

) = w

i,j

, pomijając błąd odcięcia, ostatecznie mamy:

dla i = 1,2, .. m-1

:

(46)

Powyższy wzór pozwala na określenie warunków początkowych dla j = 1, wadą tego jest
duży błąd oszacowania rzędu O(k). Alternatywą do tego jest wzór wyprowadzony jako
rozwiniecie funkcji u(x

i

,t

1

) w szereg Maclurina 2 rzędu:

dla i = 1,2, .. m-1

(47)

gdzie:

:

l

x





0

)

(

)

0

,

(

x

g

x

t

u







)

~

,

(

2

)

0

,

(

)

,

(

)

0

,

(

2

2

1

j

i

i

i

x

t

u

k

k

x

u

t

x

u

x

t

u

















)

,

0

(

~

1

t

j





)

~

,

(

2

)

0

,

(

)

0

,

(

)

,

(

2

2

1

j

i

i

i

x

t

u

k

x

t

u

k

x

u

t

x

u

















)

~

,

(

2

)

(

)

0

,

(

)

,

(

2

2

1

j

i

i

i

i

x

t

u

k

x

g

k

x

u

t

x

u













)

(

0

,

1

,

i

i

i

x

g

k

w

w









)

(

2

)

1

(

0

,

1

0

,

1

2

0

,

2

1

,

i

i

i

i

i

x

g

k

w

w

w

w



















h

k



 

110

Równania hiperboliczne - Falowe - przykład 1

Rozwiązać poniższe równanie falowe:

dla

oraz

z warunkami brzegowymi

dla

i początkowymi

oraz

dla

Obliczenia wykonać dla czasu t = 1 wykorzystując metodę MRS dla równań falowych,
przyjąć: h = 0.1 oraz k = 0.05 ,przyjąć



2

= 4

Rozwiązanie dokładne ma postać:

Tabela wyników

0

)

,

(

)

,

(

2

2

2

2

2





t

x

dx

u

d

t

x

dt

u

d



0

)

,

1

(

)

,

0

(





t

u

t

u

1

0





x

0



t

0



t

)

sin(

)

0

,

(

x

x

u





1

0





x

0

)

0

,

(







x

t

u

)

2

cos(

)

sin(

)

,

(

t

x

t

x

u







0

0,8090169944

0,8090169944

0,3

0

0,3090169944

0,5877852523

0,8090169944

0,9510565163

1,0000000000

0,9510565163

0,5877852523

0,3090169944

0

u( x

i

, 1

)

5,421



10

-20

0

0

1,1102



10

-16

1,1102



10

-16

0

0

0

| u( x

i

, 1

)

– w

i,20

|

0

0,3090169944

0,5877852523

0,8090169944

0,9510565163

1,0000000000

0,9510565163

0,5877852523

0,3090169944

0

w

i,20

k = 0.05

0,5

0,6

0,7

0,8

0,9

1,0

0,4

0,2

0,1

0,0

x

i

111

Interpretacja graficzna

112

Równania hiperboliczne – kody źródłowe
Funkcja określająca warunki początkowe równania paraboliczne

function F( X,Y : real ) : real;

begin
F := sin( pi * x );
end;

function G( X,Y : real ) : real;
begin
G := 0;
end;

Definicja wymaganych stałych i typów

const
max = 100;
type
Vector = array [ 0..max ] of real;
Macierz = array [ 0..max ] of Vector;

Funkcja „WaveEquation” – oblicza przybliżenie rozwiązania równania falowego

function WaveEquation(FX,ALPHA ,FT : real; M,N : integer; var W : Macierz):
boolean ;
var H,K,V : real; M1,N1,I,J : integer;
begin
M1 := M + 1; N1 := N + 1; Result := True;
H := FX / M; K := FT / N;
V := ALPHA * K / H;
for J := 1 to N do begin
W[0,J] := 0; W[M,J] := 0;
end;
W[0,0] := F( 0 ); W[M,0] := F ( FX );
for I := 1 to M-1 do begin
W[I,0] := F( H * ( I ) );
W[I,1] := (1-V*V)*F(H*(I))+V*V*(F((I+1)*H) + F(H*(I-1)))/2+K*G(H*(I));
end;
for J := 1 to N-1 do for I := 1 to M-1 do
W[I,J+1] := 2*(1-V*V)*W[I,J]+V*V*(W[I+1,J]+W[I-1,J])-W[I,J-1];
end;