Różniczkowanie numeryczne

(1)

Różniczkowanie numeryczne

Różniczkowaniem numerycznym nazywamy wyznaczanie przybliżonych wartości pochodnych funkcji dyskretnej jednej lub wielu zmiennych w zadanych punktach obszaru. Operację taką można wykonać dwuetapowo:

• wyprowadzając stosowne wzory różnicowe pozwalające na wyrażenie pożądanej pochodnej w określonym punkcie obszaru jako kombinacji liniowej wartości funkcji w punktach są- siednich – węzłach (tzw. schematy różnicowe),

• stosując wyznaczone wzory różnicowe do obliczania pochodnych w wybranych punktach obszaru.

Postępowanie takie jest wyjątkowo skuteczne i efektywne w przypadku, gdy węzły są rozmieszczone równomiernie, gdyż wówczas wystarczy jednokrotne wyprowadzenie schematu różnicowego na każdą niezbędną pochodną.

Metody generacji wzorów różnicowych możemy podzielić na trzy kategorie:

• przez interpolację,

• nieoznaczonych współczynników,

• ścisłe dla wielomianów możliwie najwyższego stopnia,

ponadto, w przypadku generowania wzorów różnicowych na pochodne rzędu wyższego niż pierw- szy możemy zastosować metodę składania operatorów, pamiętając że druga pochodna to pierwsza pochodna z pierwszej pochodnej, z kolei trzecia pochodna to druga pochodna z pierwszej pochodnej, lub pierwsza pochodna z drugiej pochodnej, itd.:

M

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛ ⎟

⎠

⎜ ⎞

⎝

= ⎛

⎟⎟⎠

⎜⎜ ⎞

⎝

= ⎛

⎟⎠

⎜ ⎞

⎝

= ⎛

⎟⎠

⎜ ⎞

⎝

= ⎛

dx d dx

d dx

d

dx d dx

d dx

d

2 2 2

2 3 3

2 2

. (1)

Metoda generacji wzorów różnicowych przez interpolację sprowadza się do utworzenia właściwego wielomianu interpolacyjnego (patrz rozdział poświęcony interpolacji):

( ) ( ) ∑ ( )

=

⋅

=

≅ ⁿ

i i i

n x a x

P x f

0

ϕ , (2)

a następnie jego zróżniczkowania stosowną liczbę razy:

( ) ( ) ( )

( ) ( ) ∑ ( )

∑

=

ϕ

⋅

=

≅

ϕ

⋅

=

≅

n

i i i

n

i i i

n

x a x

P x f

x a x

P x f

0

"

0 ' '

"

'

. (3)

Należy tu zwrócić uwagę, że przy takim trybie postępowania różniczkuje się wyłącznie funkcje ba- zowe, natomiast współczynniki interpolacji nie ulegają zmianie.

Metoda nieoznaczonych współczynników bazuje na rozwinięciu funkcji w szereg Taylora w otoczeniu punktu, dla którego konstruowany będzie wzór różnicowy i zażądaniu poprawności odtwo- rzenia wartości poszukiwanej pochodnej jako kombinacji liniowej wartości funkcji w węzłach ota- czających tenże punkt.

(2)

Najprostsza koncepcyjnie jest trzecia metoda, sprowadzająca się do egzekwowania wymogu, aby tworzony wzór różnicowy dawał wynik ścisły dla jednomianów możliwie wysokiego stopnia.

Poniżej przedstawimy sposób generacji wzoru różnicowego na pierwszą i drugą pochodną przy zastosowaniu każdej z wymienionych powyżej metod. Przyjmiemy przy tym, że węzły są równomier- nie rozmieszczone w odległości h jeden od drugiego, a wzór różnicowy będzie wyprowadzony dla węzła centralnego w gwieździe.

Metoda przez interpolację

Rozważamy wielomian interpolacyjny Lagrange’a przechodzący przez trzy węzły x₀ ,x₁ ,x₂:

( ) ( ) ( )

( ) ( ) ( ) ( )

( ) ( )

( ) ( ) ( ) ( ) (

₂ ⁰₀

) (

₂¹ ₁⁰

)

¹

2 2 2

1 0 1

2 0 2 0 2 1 2

0 1 0

2 1 2 1 2 0

1 2 0 2

1 0

2 2 1 0 1

2 0

1 2 0 1 0

2 1

0 2

x x x x

x x x x x x f x

x x x x

x x x x x x f x

x x x x

x x x x x x f x

x x x x

x x x f x

x x x x

x x x f x

x x x x

x x x f x

x P

−

⋅

−

⋅ +

⋅

−

⋅

⋅ −

− +

⋅

−

⋅ +

⋅

−

⋅

⋅ −

− +

⋅

−

⋅ +

⋅

−

⋅

⋅ −

=

−

⋅

−

⋅

⋅ −

− +

⋅

−

⋅

⋅ −

− +

⋅

−

⋅

⋅ −

=

, (4)

po jego zróżniczkowaniu otrzymujemy:

( ) ( ) ( ) ( ) ( ) ( ) ( )

( )

⁰

(

₀ ₁

) (

₀ ₂

)

¹

(

₁ ₀

) (

₁ ₂

)

²

(

₂ ₀

) (

₂ ₁

)

"

2

1 2 0 2

1 2 0

2 1 0 1

2 1 0

2 0 1 0

2 0 1

' 2

2 2

2

2 2

2

x x x f x

x P

x x x x

x x f x

x x x x

x x f x

x x x x

x x f x

x P

−

⋅

⋅ −

− +

⋅

⋅ −

− +

⋅

⋅ −

=

−

⋅

−

⋅ ⋅

− +

⋅

−

⋅ ⋅

− +

⋅

−

⋅ ⋅

=

, (5)

co w przypadku węzłów równoodległych prowadzi do:

( )

₂ ₀ ₂ ₁ ₂ ₂

"

2

2 1 2 1 0

20 2

2 2 0 ' 1

2

1 2

2 2

h f h f

h f x

P

h f x x f x

h x x f x

h x x x x

P

⋅ +

⋅

−

⋅

=

⋅ ⋅

−

− + ⋅

− ⋅

−

− ⋅

⋅ ⋅

−

= ⋅

. (6)

Zgodnie z zależnością (6)2 wzór różnicowy na drugą pochodną jest identyczny bez względu na po- łożenie punktu, w którym tę pochodną chcemy obliczyć, natomiast wzór (6)¹ prowadzi do dobrze znanych wzorów różnicowych dla punktów lewego, centralnego i prawego w schemacie różnico- wym (gwieździe):

( ) ( )

( )

₂ ₀ ₁ ₂

' 2 2

2 0

1 ' 2 1

2 1

0 0

' 2 0

2 3 2

2

1 2

2 3

h f h f

h f x

P x x

h f h f

x P x x

h f h f

h f x

P x x

⋅ ⋅

−

⋅ +

⋅ ⋅

−

=

⇒

=

⋅ ⋅ +

⋅ ⋅

−

=

⇒

=

⋅ ⋅

−

⋅ +

⋅ ⋅

−

=

⇒

=

. (7)

Metoda nieoznaczonych współczynników Skorzystamy z rozwinięcia funkcji w szereg Taylora:

( ) ( ) ( ) ( ) ( ) ( )

^' ² ^''

( )

( ) ( )

^{( )}¹

... !

! 2

! 1

− +

+

− ⋅ + +

− ⋅ +

= ⁿ f ⁿ a O x a ⁿ

n a a x

a f a x

f x

f ..., (8)

który pozwala wyrazić wartość funkcji w wybranym punkcie jako kombinację liniową wartości funkcji i jej kolejnych pochodnych w innym punkcie leżącym w jego sąsiedztwie.

Dążymy do wyrażenia wartości pochodnej w pożądanym punkcie jako kombinacji liniowej warto- ści funkcji w węzłach otaczających, czyli w przypadku środkowego węzła siatki trzech równomier-

(3)

nie rozmieszczonych węzłów (oznaczonych odpowiednio indeksami i−1 , i , i+1), i wzorów róż- nicowych na pierwszą i drugą pochodną otrzymamy:

1 1 1

1 '' '

1 1 1

1 '' '

1 0 0

0 1 0

+ +

−

+ +

−

⋅ β +

⋅ β

=

⋅ +

⋅

⋅ α +

⋅ α

=

⋅ +

⋅

i i i i i i i i i

f f

f

f . (9)

Rozwijając funkcję f w szereg Taylora w otoczeniu węzła i dla każdego z węzłów x_i₋₁ , x_i , x_i₊₁ (czyli we wzorze (8) przyjmując odpowiednio a= , i kolejno x_i x=x_i₋₁,x_i,x_i₊₁ , ponadto uwzględ- niając, że x_i₋₁−x_i =−h a x_i₊₁−x_i=h) oraz dokonując obcięcia rozwinięcia na trzecim wyrazie sze- regu otrzymujemy dla wzoru różnicowego na pierwszą pochodną:

1 2 ''

' 1

1 '' 2 ' 1

! 2

! 1

! 2

! 1

+ +

−

α

⋅

⋅ +

= +

α

⋅

=

α

⋅

⋅ +

⋅

−

=

i i i

i i

i

i i i

i i

h f h f

f f

h f h f

f f

, (10)

a dla wzoru różnicowego na drugą pochodną:

1 2 ''

' 1

1 2 ''

' 1

! 2

! 1

! 2

! 1

+ +

−

β

⋅

⋅ +

= +

β

⋅

=

β

⋅

⋅ +

⋅

−

=

i i i

i i

i

i i i

i i

h f h f

f f

h f h f

f f

. (11)

Po wykonaniu dodawania stronami i przegrupowaniu składników w sumie po prawej stronie mo- żemy ostatecznie zapisać w przypadku pierwszym i drugim:

( ) ( ) ( )

(

₁ ₁

)

^'

(

₁ ₁

)

^'' ²

(

₁ ₁

)

1 1 1

1

1 1 '' 2

1 1 '

1 1

1 1 1

1

! 2

! 1

! 2

! 1

+

−

− + +

− +

+

−

+

−

− + +

− +

+

−

β + β

⋅ + β

− β

⋅ + β + β + β

⋅

=

⋅ β +

⋅ β

α + α

⋅ + α

− α

⋅ + α + α + α

⋅

=

⋅ α +

⋅ α

i i i

i i i i i i i i i i i i i

i i i

i i i i i i i i i i i i i

f h f h

f f f

f

f h f h

f f f

f

. (12)

Żądając, aby wyrażenia stojące po lewej stronie wzorów (9) były odpowiednio równe wyrażeniom po prawej stronie wzorów (12) dla dowolnych wartości funkcji f i jej pochodnych, otrzymujemy dwa niewielkie układy liniowych równań algebraicznych na wyznaczenie nieznanych współczynni- ków wzorów różnicowych:

( )

⎪⎪

⎪

⎩

⎪⎪

⎪

⎨

⎧

= ⋅

= +

⋅

=

⇒

=

−

⋅

= −

= + +

+ +

−

− +

−

h h

i i

i

i i

i

i i

2 0 1

! 2

0

! 1 1

2 0 1

1 1

1 2

1 1

1

α α

α

α α

α

α α

, (13)

( )

⎪⎪

⎪

⎩

⎪⎪

⎪

⎨

⎧

=

= +

⋅

−

=

⇒

=

−

⋅

=

= + +

+ +

−

− +

−

1 2 1

1 2

1 2 1

1

1 1

! 2

0 2

! 1

0 1

h

h h

i i

i

i i

i

i i

β β

β

β β

β

β β

. (14)

(4)

Metoda wzorów ścisłych dla wielomianu

W metodzie tej żądamy, aby wyprowadzony wzór różnicowy dawał wyniki ścisłe dla jednomianów możliwie najwyższego stopnia, poczynając od 1,x,x²,...,xⁿ,... , co prowadzi do układu liniowych równań algebraicznych, z którego wyznaczamy współczynniki operatora różnicowego.

Dla wzoru różnicowego na pierwszą czy drugą pochodną zbudowanego na trzech węzłach (9) mo- żemy zażądać ścisłości wyników dla wielomianów do stopnia drugiego włącznie, co prowadzi odpowiednio do równań:

h h h

h x

x

h h

x

f f f

i i i

= ⋅ α

= α=− ⋅ α

⎪⇒

⎭

⎪⎬

⎫

=

⋅ α +

⋅ α

⇒

⋅

=

⋅ α +

⋅ α

−

⇒

=

⋅ α +

⋅ α

⇒

+

−

+

−

+

−

+

−

2 1 02 1

0 0

2 2

1 0

0 1

0 1 1

1 0

0 1

1 1

2 1 2

1 2

1 1

'' '

, (15)

dla wzoru na pierwszą pochodną i:

1 2 2 1 2

2 1 2

1 2

1 1

'' '

1 2 1

2 0

2 2

0 0

0 1

0 1 1

1 0

0 1

h h h

h h

x x

h h

x

f f f

i i i

= β

−

= β

⎪⇒

⎭

⎪⎬

⎫

=

⋅ β +

⋅ β

⇒

⋅

=

⋅ β +

⋅ β

−

⇒

=

⋅ β +

⋅ β

⇒

+

−

+

−

+

−

+

− . (16)

dla wzoru na drugą pochodną.

Wzory (15) i (16) zostały zapisane w lokalnym układzie współrzędnych, zaczepionym w węźle centralnym, czyli węźle, dla którego wyznaczany jest schemat różnicowy. W węźle tym żądamy kolejno, aby pierwsza (15) i druga (16) pochodna z:

• funkcji stałej były równe zero (równanie pierwsze),

• funkcji liniowej były odpowiednio równe jeden i zero,

• funkcji kwadratowej były odpowiednio równe zero i dwa.

Pierwsze trzy kolumny we wzorach (15) i (16) oznaczają odpowiednio: stosowny jednomian, jego pierwszą i drugą pochodną, natomiast w czwartej kolumnie po lewej stronie znaku równości obli- czono wartość wyrażenia (9) dla rozpatrywanej funkcji, a po prawej stronie umieszczono wartość pochodnej w węźle centralnym (wynikającą wprost z kolumny drugiej lub trzeciej).

Należy tu zwrócić uwagę na fakt, że istnieje ścisła zależność wiążąca liczbę węzłów w schemacie różnicowym i rząd wyznaczanej pochodnej. I tak do skonstruowania schematu różnicowego na pierwszą pochodną potrzeba minimum 2, drugą 3 itd., nie pokrywających się węzłów, czyli ogólnie wzór różnicowy na pochodną n rzędu musi zawierać co najmniej n+1 węzłów w zadaniu jednowymiarowym. Zastosowanie liczby węzłów większej od przedstawionego minimum pozwala na uzyskanie wzorów różnicowych wyższej dokładności, jednakże należy pamiętać, że liczba węzłów w gwieździe nie może być przesadnie duża ze względu na występowanie efektu Rungego (patrz rozdział o interpolacji). Pamiętać należy również o tym, że dokładność schematu różnicowego jest najwyższa w pobliżu środka ciężkości zbioru węzłów, na których jest on tworzony, i jest to nieza- leżne od rzędu pochodnej. Warto również zwrócić uwagę na fakt, że wszystkie trzy metody przedstawione powyżej w przypadku takiego samego rozmieszczenia węzłów w schemacie różnicowym

(5)

prowadzą do identycznych wartości współczynników operatora różnicowego. Nie jest to wynik przypadkowy, gdyż, jak łatwo zauważyć, wszystkie metody mają swoje podstawy w interpolacji wielomianowej.

Przedstawione powyżej metody można uogólnić na przypadek dowolnej liczby węzłów w schemacie różnicowym i pochodnej dowolnego rzędu, co pokażemy poniżej na przykładzie metody nieoznaczonych współczynników, należy jednak pamiętać o tym, że każda kolejna pochodna jest obli- czana z większym błędem.

W celu wyprowadzenia wzorów ogólnych w zadaniu jednowymiarowym rozważmy zbiór węzłów ponumerowanych ,...,i−2,i−1,i,i+1,i+2 , z których węzeł i jest węzłem centralnym. Poszuku- jemy współczynników wzoru różnicowego na n – tą pochodną w tym węźle:

( )ⁿ =...+ _i₋₂⋅ _i₋₂+ _i₋₁⋅ _i₋₁+ _i⋅ _i+ _i₊₁⋅ _i₊₁+ _i₊₂⋅ _i₊₂+...

i f f f f f

f α α α α α . (17)

Postępując podobnie jak poprzednio, czyli rozwijając funkcję f w szereg Taylora (8) w otoczeniu węzła i dla każdego węzła sąsiadującego, i uwzględniając, że:

,...h_k =x_k−x_i , k=...,i−2,i−1,i,i+1,i+2 , (18) otrzymujemy:

( ) ( )

( )

M M

+

α

⋅

⋅ + +

⋅ +

=

α

⋅

⋅ + +

⋅ +

=

α

⋅

=

α

⋅

⋅ + +

⋅ +

=

α

⋅

⋅ + +

⋅ +

=

+ + +

− −

−

− −

−

− −

2 ''' 2

3 '' 2 2 ' 2 2 2

1 1 3 '''

'' 1 2 ' 1 1 1

1 ''' 1

3 '' 1 2 ' 1 1 1

2 2 3 '''

'' 2 2 ' 2 2 2

... !

! 3

! 2

! 1

... !

! 3

! 2

! 1

... !

! 3

! 2

! 1

... !

! 3

! 2

! 1

i n i n i i

i i i i i i i

i n i n i i

i i i i

i i i

i i

i

i n i n i i

i i i i i i i

i n i n i i

i i i i

i i i

n f f h

f h f h

f h f

n f f h

f h f h

f h f

f f

n f f h

f h f h

f h f

n f f h

f h f h

f h f

. (19)

Po wykonaniu dodawania stronami i przegrupowaniu składników w sumie po prawej stronie znaku równości można zapisać:

( )

(

^... ^...

)

...

'' '

...

''

...

'

...

1

2 2 1! 1 1 1! 1 1 1! 2 2 1!

2 3

2

! 31 1 3

1

! 31 1 3

1

! 31 2 3

2

! 31

2 2

2

! 21 1 2

1

! 21 1 2

1

! 21 2 2

2

! 21

2 2 1 1 1 1 2 2

2 1 1

2

+ α

⋅

⋅ + α

⋅

⋅ + α

⋅

⋅ + α

⋅

⋅ +

+ α

⋅

⋅ + α

⋅

⋅ + α

⋅

⋅ + α

⋅

⋅ +

+ α

⋅

⋅ + α

⋅

⋅ + α

⋅

⋅ + α

⋅

⋅ +

+ α

⋅ + α

⋅ +

+ α + α + α + α + α +

⋅

=

⋅

+ + +

+

−

+ + +

+

−

+ + +

+

−

+ + + +

−

+ +

−

i n i n i n i n i n i n i n i n n

i

i i i

i i

i i i

i i

i i i i i i i i i

i i i i i i

n i

h h

f

h h

f

h h

f

h h

f f f

M M

, (20)

co ostatecznie prowadzi do układu liniowych równań algebraicznych opisanego macierzą Vander- monde’a:

(6)

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢

⎣

⎡

=

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢

⎣

⎡

⋅

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

⋅

+ +

−

+ +

−

+ +

−

+ +

−

+ +

−

1 0 0 0 0

0 0 0 0 0

1 1

1 1 1

2 1 1 2

2 1! 1 1! 1 1! 2 1!

3 2

! 31 3

1

! 31 3

1

! 31 3

2

! 31

2 1

! 21 2

1

! 21 2

1

! 21 2

2

! 21

2 1

1 2

M M

L L

L L L

L L

L

L L

i i i i i

n i n n i n n i n n i n

i i

i i i

i

h h

h h h

h

h h

α α α α α

, (21)

z którego wyznaczamy nieznane współczynniki:

,...α_k , k=...,i−2,i−1,i,i+1,i+2 . (22) Dla przedstawienia związku pomiędzy liczbą węzłów gwiazdy, odległościami pomiędzy nimi i ja- kością uzyskanego numerycznego oszacowania wartości pochodnej dokonamy obliczenia pierwszej pochodnej dla skrajnego prawego węzła regularnej siatki węzłów, używając kolejno wzorów różni- cowych zbudowanych na od 2 (minimalna możliwa liczba węzłów ze względu na rząd pochodnej) do 5 węzłów. Obliczenia przeprowadzimy w punkcie x=1,0 dla funkcji y=⁵ x , której pochodna wyznaczona analitycznie jest równa

5 4

1 5 ' 1

y = ⋅ x , a jej wartość dla x=1,0 wynosi y'=0,2 .

Współczynniki operatora różnicowego na pierwszą pochodną dla poszczególnych wzorów różni- cowych wyznaczamy z układu równań (21), uwzględniając w jego zapisie liczbę węzłów gwiazdy i ich względną lokalizację oraz położenie węzła, w którym generowany jest operator. Poniżej, dla przykładu zapiszemy ten układ równań i jego rozwiązanie w przypadku najbardziej złożonym – gwiazdy pięciowęzłowej:

( ) ( ) ( ) ( )

h h

h

h h

i i i i i

⋅ ⋅

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

−

=

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

⇒

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

=

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

⋅

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

−

⋅

−

⋅

−

⋅

−

⋅

−

⋅

−

⋅

−

⋅

−

⋅

−

⋅

−

⋅

−

⋅

−

⋅

−

⋅

−

⋅

−

⋅

−

12 1

25 48 36

16 3

0 0 0 1 0

0 2

3 4

0 2

3 4

0 2

3 4

0 2

3 4

1 1 1

1 1

1 2 3 4

4

! 41 4

! 41

3

! 31 3

! 31

2

! 21 2

! 21

α α α α α

. (23)

Wyniki uzyskane dla tego i pozostałych, prostszych, przypadków zebrano w tablicy 1, która zawiera współczynniki wzoru różnicowego na pierwszą pochodną w skrajnym prawym węźle siatki (oznaczonym jako i ), zbudowanego na 2, 3, 4 i pięciu węzłach n.

Tablica 1. Współczynniki wzoru różnicowego.

n α_i₋₄ α_i₋₃ α_i₋₂ α_i₋₁ α_i

2 -1 1 × 1h

3 1 -4 3 ×¹₂⋅h

4 -2 9 -18 11 ×¹₆⋅h

5 3 -16 36 -48 25 ×¹₁₂⋅h

Współczynniki zebrane w tablicy 1 zostały wykorzystane do obliczenia przybliżonej wartości pierwszej pochodnej funkcji y=⁵ x w punkcie x=1,0 .

(7)

Tablica 2. Przybliżona wartość pierwszej pochodnej funkcji y=⁵ x w punkcie x=1,0 . Odległość między węzłami h i błąd względny procentowy pochodnej b n h=10,0 _b h=1,0 b h=0,1 _b h=0,01 _b 2 0,255185 27,59 1,000000 400,00 0,208516 4,2582 0,200805 0,4024 3 0,370270 85,13 1,000000 400,00 0,198795 0,6023 0,199999 0,0049 4 0,443957 121,98 0,716233 258,11 0,200337 0,1684 0,200000 0,0001 5 0,497927 148,96 0,303498 51,75 0,199843 0,0786 0,200000 0,0000

Jak wynika z tablicy 2, przy dostatecznie małych odległościach pomiędzy węzłami zwiększanie liczby węzłów w gwieździe zdecydowanie poprawia jakość uzyskanego oszacowania, jednak w przypadku gdy odległości pomiędzy węzłami są duże, może zdarzyć się sytuacja, kiedy oszacowanie uzyskane przy zastosowaniu wzoru niższego rzędu będzie lepsze niż oszacowanie uzyskane przy zastosowaniu wzoru zawierającego więcej węzłów. Niestety nie jest możliwe sformułowanie ogólnego wzoru na tę zależność.

Każda z metod przedstawionych wyżej może być również uogólniona na zadania dwu i trój wymia- rowe. Główną trudnością w tym przypadku jest zagwarantowanie odpowiedniego (na płaszczyźnie bądź w przestrzeni) rozmieszczenia węzłów, zapewniającego dostateczną stabilność wyprowadza- nych wzorów. Oczywiście problem ten nie istnieje w przypadku regularnych siatek węzłów. Poni- żej, dla ilustracji, przedstawimy wyprowadzenie wzoru różnicowego na operator Laplace’a w zadaniu dwuwymiarowym dla regularnej siatki węzłów. Zastosujemy metodę nieoznaczonych współ- czynników.

W zadaniu dwuwymiarowym wzór Taylora przedstawia się następująco:

( ) ( )

_{( )}

( ) ( )

( ) ( ) ( )

( ) ( ) ( ) ( )

L

+

⋅ +

⋅ ⋅ + ⋅

⋅

+

⋅ +

⋅ ⋅ + ⋅

⋅

+

⋅ +

⋅

+

=

= + +

''' , , ''' 3

, , ''' 2

, , 3

'' , , '' 2

, , ''

, , 2

' , , '

, ,

,

! 3

! 3 3

! 3

! 2

! 2 2

! 2

! 1

! 1 ,

,

yyy j i xyy

j i xxy

j i xxx

j i

yy j i xy

j i xx

j i

y j i x

j i

k f k f

f h k f h

h

k f k f

f h h

k f h f

f k

h f k y h x f

(24)

gdzie indeksy

( )

i,j oznaczają numer wiersza i kolumny regularnej prostokątnej siatki węzłów, na przecięciu których znajduje się rozpatrywany węzeł centralny.

Poszukujemy wzoru różnicowego na operator Laplace’a, czyli:

( ) ( ) ( ) ( ) ( )^'', ,

'' , , '

, , '

, , 2 ,

2 2 2

1 1

0 0

0 f_i _j f_i_j _x f_i _j _y f_i_j _xx f_i _j _yy y

x = ⋅ + ⋅ + ⋅ + ⋅ + ⋅

∂ + ∂

∂

∂ . (25)

Wobec tego zapisujemy rozwinięcie funkcji f w szereg Taylora w węźle centralnym gwiazdy i czterech najbliższych mu węzłach sąsiadujących, opatrzonych odpowiednio indeksami

(

i−1, j

)

,

(

i+1,j

)

,

(

i, j−1

)

oraz

(

i, j+1

)

:

(8)

( ) ( )

( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( ) ( ) 4

''' , ,

! 3 ''

, ,

! 2 '

, ,

! 1 , 1 ,

3 '''

, ,

! 3 ''

, ,

! 2 '

, ,

! 1 , 1 ,

2 '''

, ,

! 3 ''

, ,

! 2 '

, ,

! 1 , ,

1

1 '''

, ,

! 3 ''

, ,

! 2 '

, ,

! 1 , ,

1

0 ,

,

3 2

α α α α α

⋅ +

= +

⋅ +

⋅

−

⋅ +

⋅

−

=

⋅ +

=

⋅ +

⋅

−

⋅ +

⋅

−

=

⋅

=

+

− +

−

K K K K

yyy j k i yy j k i y j k i j i j

i

yyy j k i yy j k i y j k i j i j

i

xxx j h i xx j h i x j h i j i j i

j i j

i

f f

f

f f

f

f f

f

f f

f

f f

. (26)

We wzorach (26) uwzględniono fakt, że wybrane współczynniki rozwinięć zapisanych w wierszach 2 – 4 są równe zeru ponieważ h=0 lub k=0 . Dalej, podobnie jak poprzednio, dokonujemy ob- cięcia we wszystkich rozwinięciach po drugiej pochodnej, dodajemy poszczególne równania stronami, przegrupowujemy prawą stronę uzyskanego w ten sposób równania tak, aby przed nawiasy wyciągnąć funkcję i jej poszczególne pochodne i porównujemy ją z prawą stroną równania (25), co, po przyjęciu dodatkowo h=k =d (siatka regularna kwadratowa), prowadzi do:

( ) ( ) ( ) ( ) ( ) ( )

( )

_{( )}

( )

( )^''_, _, ²

(

₁ ₂

)

_{( )}^''_, _, ²

(

₃ ₄

)

3 4 '

, , 1 2 '

, ,

4 3 2 1 0 , ''

, , ''

, , '

, , ,

2 2

1 1

0 0

0

α α α

α

α α α

α

α α α α α

+

⋅

⋅ +

+

⋅

+

−

⋅

⋅ +

−

⋅

+ + + + +

⋅

=

⋅ +

⋅

f d f d

d f d

f f f

f f

yy j i xx

j i

y j i x

j i

j i yy

j i xx j i y j i x

j i j

i

, (27)

i ostatecznie układu liniowych równań algebraicznych na wyznaczenie niewiadomych współczyn- ników α , gdyż żądamy aby prawa i lewa strona równości (27) były sobie równe dla dowolnych wartości funkcji f i jej pochodnych występujących w tym wzorze.

( )

⎪⎪

⎩

⎪⎪

⎨

⎧

×

=

= +

⋅

=

= +

⋅

⇒ =

=

−

⋅

=

−

⋅

−

=

= + + + +

2

4 4

3 2

1

2 4 3 2

1 1

2

0 4

3 2 1 0

1

1 2 1

1 0

4 0

d d

d d d

α α

α

α α

α

α α

α

α α

α

α α

α α α α

. (28)

Wyznaczymy teraz współczynniki tego samego operatora różnicowego posługując się metodą skła- dania operatorów. W tym celu wystarczy zauważyć, że operator Laplace’a to suma drugich pochodnych w kierunku x i y , a więc do jej obliczenia wystarczy zastosować dwukrotnie na przy- kład wyprowadzony powyżej (14) wzór różnicowy na drugą pochodną w zadaniu jednowymiarowym, sumując współczynniki w pokrywających się węzłach. W rezultacie, dla siatki regularnej kwadratowej, otrzymamy wynik identyczny jak uzyskany powyżej (28), gdyż:

( ) ( ) ( ) ( )

( ) 2 (, 1) 2 ( ), 2 (, 1)

'' , 2 ,

2

, 2 1 2 ,

, 2 1 ''

, 2 ,

2

1 2

1

1 2

1

+

−

+

−

⋅ +

⋅

≅

∂ =

∂

⋅ +

⋅

≅

∂ =

∂

j i j

i j

i yy

j i

j i j

i j

i xx

j i

d f d f

d f y f

d f d f

d f

x f , (29)

przy oznaczeniach takich samych jak użyte we wzorze (27), a po zsumowaniu:

( ) ( ) 2 ( 1, ) 2 (, 1) 2 ( ), 2 ( 1, ) 2 (, 1)

'' , , ''

, 2 ,

2 2

2 1 1 4 1 1

+ +

−

− + ⋅ − ⋅ + ⋅ + ⋅

⋅

≅ +

∂ = + ∂

∂

j i j

i j

i yy

j i xx j

i f

f d f d

f d y f

x . (30)