Metody numeryczne rozwiązywania równań

różniczkowych cząstkowych

Metoda różnic skończonych (siatek)

Uwagi ogólne

Dane równanie różniczkowe cząstkowe opisane operatorem L:

f

Lu

w obszarze  i warunki brzegowe:

1,2,...p

i

dla

u

l

i

x

i

i















X

k

W metodach różnicowych poszukuje się tablicy wartości
przybliżonych u

h

rozwiązania dokładnego u na zbiorze

izolowanych punktów X

k

 (k=1,2,...,N

h

)

zwanym siatką.

Punkty X

k

są nazywane węzłami siatki.

węzeł pomocniczy

węzeł podstawowy

Równania służące do wyznaczania wartości przybliżonych

nazywamy równaniami różnicowymi.

x

y

h

x

h

y

h=(h

x

,h

y

)

Parametr h
charakteryzuje
siatkę 

h

Dla równania różniczkowego:

f

Lu

w obszarze  z warunkami brzegowymi:

1,2,...p

i

dla

u

l

i

x

i

i











otrzymujemy jego odpowiednik różnicowy:

ih

h

ih

h

h

h

u

r

f

u

R







Zakładając, że problem opisany równaniem różniczkowym ma
jednoznaczne rozwiązanie, to równania różnicowe będą jego
odpowiednikiem jeżeli są spełnione następujące warunki:

ih

h

ih

h

h

h

u

r

f

u

R







1. Układ równań różnicowych posiada jednoznaczne rozwiązanie:

 





h

h

k

h

k

h

X

dla

X

u





dla każdego dopuszczalnego h.

2. Zbieżność do rozwiązania dokładnego u.

Oznacza to, że rozwiązanie u

h

powinno przy h 0 dążyć do

rozwiązania dokładnego u.

Dla określenia zbieżności jest koniecznym wprowadzenie

odpowiednich przestrzeni funkcyjnych i norm w nich.

Wprowadzamy przestrzeń funkcyjną U z normą || ||

U

, do której

należy rozwiązanie dokładne u.

oraz przestrzeń N

h

- wymiarową U

h

z normą || ||

Uh

,

której elementami są układy N

h

liczb i do której

należy rozwiązanie u

h

.

Normy || ||

U

i || ||

Uh

winny być zgodne, tzn. ponieważ funkcja u(x)

jest określona w węzłach podstawowych X

k

siatki, to mówimy,

że normy są zgodne jeżeli zachodzi:

 

 

U

Uh

0

h

x

u

x

u

lim





dla każdego uU.

Przykłady norm zgodnych:

 











d

f

f

2

L

2











0

h

h

2

j

,

i

2

ij

2

L

f

h

f

- zbiór węzłów wewnętrznych.

0

h



 

















d

,

u

u

u

i

2

i

2

H

1

2





























































h

ij

1

h

2

j

,

i

2

ij

1

j

,

i

2

ij

j

,

1

i

2

2

H

h

h

u

u

h

u

u

u

h

u

Wielkość:

 

Uh

h

h

x

u

u 





nazywamy błędem rozwiązania przybliżonego u

h

.

U

h

dąży do rozwiązania dokładnego u(x), jeżeli

 

0

x

u

u

lim

Uh

h

0

h







Jeżeli można znaleźć taką funkcję (h), że

 

 

h

x

u

u

Uh

h







to mówimy, że zostało znalezione oszacowanie błędu.

3. Stabilność

ih

h

ih

h

h

h

u

r

f

u

R







Różnicowe zagadnienie brzegowe jest stabilne, jeżeli istnieją
takie liczby >0 i h

0

>0, że dla dowolnych h<h

0

i f

h

F

h

, takich,

że ||f

h

||

Fh

< zadanie brzegowe:

ih

ih

h

ih

h

h

h

h

z

r

f

f

z

R

















posiada jedno jednoznaczne rozwiązanie i zachodzi:

Fh

h

Uh

h

h

f

C

u

z







gdzie C stała niezależna od h.

Twierdzenie wiążące stabilność i zbieżność:

Jeżeli zadanie różnicowe

ih

h

ih

h

h

h

u

r

f

u

R







jest przybliżeniem różniczkowego problemu brzegowego:

f

Lu

1,2,...p

i

dla

u

l

i

x

i

i











i rozwiązanie u

h

jest stabilne wtedy zachodzi

 

0

x

u

u

lim

Uh

h

0

h







i rząd zbieżności w funkcji h jest taki sam jak rząd aproksymacji.

Zastępowanie pochodnych ilorazami różnicowymi

na siatce prostokątnej

i-1

i

i+1

k-1

k

k+1

24

h

x

u

6

h

x

u

2

h

x

u

h

x

u

u

u

4

i

ik

4

4

3

i

ik

3

3

2

i

ik

2

2

i

ik

ik

k

,

1

i





























24

h

x

u

6

h

x

u

2

h

x

u

h

x

u

u

u

4

i

ik

4

4

3

i

ik

3

3

2

i

ik

2

2

i

ik

ik

k

,

1

i





























h

i

h

k

 
 

2

k

2

k

1

k

,

i

1

k

,

i

2

i

2

i

k

,

1

i

k

,

1

i

h

O

h

2

u

u

y

u

h

O

h

2

u

u

x

u





























i-1

i

i+1

k-1

k

k+1

h

i

h

k

Druga pochodna dodając

stronami:

24

h

x

u

6

h

x

u

2

h

x

u

h

x

u

u

u

4

i

ik

4

4

3

i

ik

3

3

2

i

ik

2

2

i

ik

ik

k

,

1

i





























24

h

x

u

6

h

x

u

2

h

x

u

h

x

u

u

u

4

i

ik

4

4

3

i

ik

3

3

2

i

ik

2

2

i

ik

ik

k

,

1

i





























 

 

2

k

2

k

1

k

,

i

ik

1

k

,

i

2

2

2

i

2

i

k

,

1

i

ik

k

,

1

i

2

2

h

O

h

u

u

2

u

y

u

h

O

h

u

u

2

u

x

u

































i-1

i

i+1

k-1

k

k+1

h

i

h

k

Dla pochodnej mieszanej



























































i

3

k

k

3

i

k

i

1

k

,

1

i

1

k

,

1

i

1

k

,

1

i

1

k

,

1

i

2

h

h

,

h

h

max

O

h

h

4

u

u

u

u

y

x

u

Konstrukcja warunków brzegowych na siatce

i

k

 

k

,

i

h



 





i



k







i

ik

x











k

ik

x







1. Przeniesienie wartości:

Przyjmujemy:

 











































i

k

k

ik

k

i

i

ik

h

dla

x

dla

x

k

,

i

r(i,k,)

lub

 

 























,

k

,

i

r

min

x

h

x

k

,

i

2. Interpolacja liniowa dla warunku brzegowego Dirichleta

i,k

i-1,k

i+1,k

h



Zakładając liniowy rozkład rozwiązania między sąsiednimi
węzłami mamy:

)

x

(

u

h

x

x

u

h

x

x

u

1

i

ik

i

k

,

1

i













i dla x=x

i

+ mamy:



 

i

ik

k

,

1

i

h

1

u

h

u





























 

i





3. Interpolacja liniowa dla warunku brzegowego Neumanna.

n

i,k

i-1,k

i,k-1

h

k

h

i



 

k

,

i

sin

h

u

u

cos

h

u

u

k

1

k

,

i

ik

i

k

,

1

i

ik

















































y

,

x

y

,

x

u

y

,

x

d

n

u

y

,

x

c

Równania eliptyczne

Dla uproszczenia rozważań będą analizowane przypadki

dwuwymiarowe x=(x,y)

















y

,

x

f

u

y

,

x

g

y

u

y

,

x

b

y

x

u

y

,

x

a

x







































Warunki brzegowe:

1

2

3

4

5

6

7

8

9 x(i)

1

2

3

4

5

6

7

8

9

10 y(k)

h

y

h

x

(0,0)

j

j+1

j-1

l

p

















y

,

x

f

u

y

,

x

g

y

u

y

,

x

b

y

x

u

y

,

x

a

x







































Dla węzłów wewnętrznych będzie:

j

j+1

j-1

l

p

j

j

j

y

l

p

y

l

p

2

y

p

j

l

j

x

1

j

1

j

x

1

j

1

j

2

x

1

j

j

1

j

j

f

u

g

h

2

u

u

h

2

b

b

h

u

u

2

u

b

h

2

u

u

h

2

a

a

h

u

u

2

u

a











































lub w formie macierzowej:

f

u

A

w



h

y

h

x

(0,0)

j

j+1

j-1

l

p

styczna

normalna































y

,

x

y

,

x

u

y

,

x

d

n

u

y

,

x

c

p

Przyjmując:

(xk

p

,yk

p

)













p

p

p

p

p

p

p

p

p

yk

,

xk

d

d

yk

,

xk

c

c

yk

,

xk











p-1

m

p



otrzymujemy:

p

p

p

p

y

p

m

p

x

1

p

p

p

u

d

sin

h

u

u

cos

h

u

u

c































Uwaga dotycząca błędu obliczeń.
Generalnie jeżeli węzły nie leżą na krzywej brzegowej
i liczymy metodą przeniesienia wartości, to dokładność
obliczeń jest rzędu h. Jeżeli węzły na krzywej brzegowej
bądź wyliczamy wartości funkcji brzegowej interpolując
liniowo dokładność wzrasta do h

2

.

W zagadnieniu Dirichleta oprócz trudności z wyznaczeniem

wartości brzegowych nie ma innych problemów i otrzymany

układ równań algebraicznych najczęściej można rozwiązać

bez kłopotów.

Sytuacja może się komplikować przy zagadnieniu Neumanna.

Siatki praktycznie nie stosowane w zagadnieniach

eliptycznych liniowych i nieliniowych.

Równania opisujące ewolucję układu w czasie

Równania paraboliczne

Dane jednowymiarowe równanie przewodnictwa:

 

 





 

 

 

 

 

 

 





T

0,

t

dla

t

f

t

1,

u

t

f

t

,

0

u

0,1

x

dla

x

0

,

x

u

T

0,

t

0,1

x

dla

t

,

x

f

,

u

,

u

2

1

0

xx

t























x

t

0

1

h



i-1 i i+1

k+1
k
k-1

x

t

0

N

h



i-1 i i+1

k+1
k
k-1

Oznaczamy operator różnicy II rzędu:

2

k

1

i

k

i

k

1

i

k

i

h

u

u

2

u

u













i wprowadzamy schematy różnicowe z wagą :









k

i

k

i

1

k

i

k

i

1

k

i

u

1

u

u

u





























N

,

0

i

i





K

,

0

k

















5

.

0

k

,

x

f

i

k

i

 

 





 

 

 

 

 

 

 





T

0,

t

dla

t

f

t

1,

u

t

f

t

,

0

u

0,1

x

dla

x

0

,

x

u

T

0,

t

0,1

x

dla

t

,

x

f

,

u

,

u

2

1

0

xx

t































k

i

k

i

1

k

i

k

i

1

k

i

u

1

u

u

u

























Problem brzegowy jest aproksymowany przez

 

i

0

i

ih

u









dla i=1,2,...N

 

 

k

2

k

N

k

1

k

0

2

1

f

k

f

u

f

k

f

u













dla k=1,2,...,K.

Warunki zgodności

k+1
k

i-1 i i+1

Schemat sześciowęzłowy









k

i

k

i

1

k

i

k

i

1

k

i

u

1

u

u

u

























Jeżeli =0 schemat jest

nazywany jawnym lub
explicite

k+1
k

i-1 i i+1

Jeżeli 0 schemat jest

nazywany niejawnym lub
implicite

k+1
k

i-1 i i+1

k+1
k

i-1 i i+1

Wartości w warstwie k+1
otrzymujemy rozwiązując
układ równań:





k

i

k

i

k

i

1

k

i

1

k

i

u

1

u

1

u

u

1





























Czysto niejawny schemat:

k+1
k

i

1





k+1
k

i-1 i i+1

Schemat Cranka - Nicholsona:

5

.

0





Oszacowanie dokładności aproksymacji.

Rozwiązanie dokładne zagadnienia brzegowego

 

 





 

 

 

 

 

 

 





T

0,

t

dla

t

f

t

1,

u

t

f

t

,

0

u

0,1

x

dla

x

0

,

x

u

T

0,

t

0,1

x

dla

t

,

x

f

,

u

,

u

2

1

0

xx

t























jest u(x,t) i jego wartość w węzłach (x

i

,t

k

) siatki będzie oznaczana

u(i,k).

Rozwiązanie zagadnienia brzegowego sformułowanego dyskretnie









k

2

k

N

k

1

k

0

i

0

i

k

i

k

i

1

k

i

k

i

1

k

i

f

u

f

u

u

u

1

u

u

u



































jest

k

i

u

Błąd aproksymacji

k

i

z

jest

)

k

,

i

(

u

u

z

k

i

k

i





Dla oceny błędu

k

i

z

w kroku k-tym wprowadza się normę, np.:





k

i

N

,

0

i

z

max

z





lub

 









1

N

1

i

2

k

i

z

h

z

Z

)

k

,

i

(

u

u

z

k

i

k

i





wynika, że

)

k

,

i

(

u

z

u

k

i

k

i





i podstawiając do
w miejsce

otrzymujemy

równoważne

zadanie różnicowe

dla









k

2

k

N

k

1

k

0

i

0

i

k

i

k

i

1

k

i

k

i

1

k

i

f

u

f

u

u

u

1

u

u

u



































k

i

u

k

i

z









0

z

0

z

0

z

z

1

z

z

z

k

N

k

0

0

i

k

i

k

i

1

k

i

k

i

1

k

i









































0

z

0

z

0

z

z

1

z

z

z

k

N

k

0

0

i

k

i

k

i

1

k

i

k

i

1

k

i

































gdzie



 

  







  

k

i

k

i

k

,

i

u

1

k

,

i

u

k

,

i

u

1

1

k

,

i

u





























błąd schematu różnicowego
w stosunku do rozwiązania
dokładnego u(x,t).









k

2

k

N

k

1

k

0

i

0

i

k

i

k

i

1

k

i

k

i

1

k

i

f

u

f

u

u

u

1

u

u

u



































Mówimy, że









k

2

k

N

k

1

k

0

i

0

i

k

i

k

i

1

k

i

k

i

1

k

i

f

u

f

u

u

u

1

u

u

u



































przybliża rozwiązanie problemu brzegowego

 

 





 

 

 

 

 

 

 





T

0,

t

dla

t

f

t

1,

u

t

f

t

,

0

u

0,1

x

dla

x

0

,

x

u

T

0,

t

0,1

x

dla

t

,

x

f

,

u

,

u

2

1

0

xx

t























z dokładnością rzędu (m,n) lub O(h

m

+

n

), jeżeli



 

  







  

k

i

k

i

k

,

i

u

1

k

,

i

u

k

,

i

u

1

1

k

,

i

u





























spełnia nierówność:





n

m

h

M









gdzie M - stała.

Dla uproszczenia zapisu wprowadzamy:



  











u

k

,

i

u

1

k

,

i

u

i mamy:







  







  







  

































u

5

.

0

k

,

i

u

1

k

,

i

u

5

.

0

k

,

i

u

1

k

,

i

u

5

.

0

k

,

i

u

1

k

,

i

u

5

.

0

1

k

,

i

u

 



  

















u

5

.

0

k

,

i

u

1

k

,

i

u

5

.

0

k

,

i

u



 

  



  





 

 



  

























































u

5

.

0

k

,

i

u

1

k

,

i

u

5

.

0

k

,

i

u

u

k

,

i

u

k

,

i

u

1

k

,

i

u

k

,

i

u

1

1

k

,

i

u

Uwzględniając powyższe równości i podstawiając do



 

  







  

k

i

k

i

k

,

i

u

1

k

,

i

u

k

,

i

u

1

1

k

,

i

u





























mamy



  









k

i

k

i

u

5

.

0

k

,

i

u

1

k

,

i

u

5

.

0



























Rozwijając funkcje u(n,p) w szereg Taylora w otoczeniu

punktu: x

i

,t

k

+0.5  oraz wprowadzając oznaczenie:





u

5

.

0

t

,

x

u

k

i







Będziemy mieli:

 

 





 

 

 









 

 

2

t

3

tt

2

3

tt

2

t

3

tt

2

t

4

xxxx

2

xx

O

,

u

u

O

,

u

8

u

)

k

,

i

(

u

1

k

,

i

u

5

.

0

O

,

u

8

,

u

5

.

0

u

k

,

i

u

O

,

u

8

,

u

5

.

0

u

1

k

,

i

u

h

O

,

u

12

h

,

u

k

,

i

u





























































Uwzględniając powyższe zależności można



  









k

i

k

i

u

5

.

0

k

,

i

u

1

k

,

i

u

5

.

0



























zapisać w postaci:













4

2

xxxx

2

txx

k

i

t

xx

k

i

h

O

,

u

12

h

,

u

5

.

0

,

u

,

u



























Ale

0

f

,

u

,

u

t

xx







gdyż u jest rozwiązaniem dokładnym

 

 





 

 

 

 

 

 

 





T

0,

t

dla

t

f

t

1,

u

t

f

t

,

0

u

0,1

x

dla

x

0

,

x

u

T

0,

t

0,1

x

dla

t

,

x

f

,

u

,

u

2

1

0

xx

t























i w każdym punkcie obszaru spełnia równanie paraboliczne.

Uwzględniając, że

xx

txx

xxxx

,

f

,

u

,

u





mamy









2

4

xx

2

txx

2

k

i

k

i

h

O

,

f

12

h

,

u

12

h

5

.

0

f















































2

4

xx

2

txx

2

k

i

k

i

h

O

,

f

12

h

,

u

12

h

5

.

0

f







































Jeżeli wyrażenie w nawiasie kwadratowym jest równe zeru, tzn.:





























12

h

5

.

0

0

12

h

5

.

0

2

2

oraz

xx

2

k

i

,

f

12

h

f 





W obliczeniach numerycznych wygodniej przyjąć:





5

.

0

k

1

i

5

.

0

k

1

i

5

.

0

k

i

k

i

f

f

12

1

f

6

5



















to schemat ma
dokładność O(h

4

+

2

)









k

2

k

N

k

1

k

0

i

0

i

k

i

k

i

1

k

i

k

i

1

k

i

f

u

f

u

u

u

1

u

u

u







































Schemat:

z





5

.

0

k

1

i

5

.

0

k

1

i

5

.

0

k

i

k

i

f

f

12

1

f

6

5



















jest schematem o podwyższonej dokładności wynoszącej:





2

4

h

O





Jeżeli  =0.5 jak w schemacie Cranka-Nicholsona, to









2

4

xx

2

txx

2

k

i

k

i

h

O

,

f

12

h

,

u

12

h

5

.

0

f















































2

4

xx

2

txx

2

k

i

k

i

h

O

,

f

12

h

,

u

12

h

5

.

0

f







































i dla =0.5 mamy:





2

2

k

i

k

i

h

O

f















Dla zachowania oceny zbieżności O(h

2

+

2

) należy przyjąć:

5

.

0

k

i

k

i

f







lub





k

i

1

k

i

k

i

f

f

5

.

0









Jeżeli 0.5 i 

*

, to dokładność obliczeń jest rzędu O(h

2

+).

Stabilność

 





 

 

 

 





T

0,

t

dla

0

t

1,

u

0

t

,

0

u

0,1

x

dla

0

0

,

x

u

T

0,

t

0,1

x

dla

,

u

,

u

xx

t



























0

u

0

u

0

u

u

1

u

u

u

k

N

k

0

0

i

k

i

1

k

i

k

i

1

k

i





























Zbadamy zachowanie się schematu:

jawnego tj.  =0





0

u

0

u

0

u

u

u

2

u

h

u

u

k

N

k

0

0

i

k

1

i

k

i

k

1

i

2

k

i

1

k

i

























Schemat jest

Przyjmijmy:

1

h

2





i załóżmy, że warunek początkowy w

punkcie i-tym jest dany z błędem  . Badamy jak przenosi się błąd
na siatce.

x

t

0

N

0

0

u

0

u

i

n

i

n

0

u

u

u

u

u

k

N

k

0

0

n

k

1

i

k

i

k

1

i

1

k

i















































 2



3



 2







 3



6



 7



6



 3







 4



10



 16



19



 16



10



 4







 5



15



 30



45



 51



45



 30



15



 5

0

0



 6



21



 50



90



 126



141



 126



90



 50



20 0

Schemat jawny
z =0 i

1

h

2





x

t

0

N

0

0

u

0

u

i

n

i

n

0

u

u

1

.

0

u

8

.

0

u

1

.

0

u

k

N

k

0

0

n

k

1

i

k

i

k

1

i

1

k

i





































1

.

0



8

.

0



1

.

0



01

.

0



16

.

0



66

.

0



16

.

0



01

.

0



01

.

0



02

.

0



2

.

0



56

.

0



2

.

0



02

.

0



01

.

0

0



01

.

0



04

.

0



22

.

0



49

.

0



22

.

0



04

.

0



01

.

0

0

0

0



01

.

0



06

.

0



23

.

0



44

.

0



23

.

0



06

.

0



01

.

0

0

0

0

0

0



01

.

0



07

.

0



23

.

0



4

.

0



07

.

0



01

.

0

0

0

Schemat jawny
z =0 i

1

.

0

h

2





Obliczenia z dokładnością do 2 miejsc



23

.

0

Analiza stabilności metodą spektralną

Niech rozwiązanie jednorodnego zagadnienia różnicowego:









0

u

0

u

u

u

1

u

u

u

k

N

k

0

n

0

n

k

n

1

k

n

k

n

1

k

n































ma postać











in

exp

u

k

k

n

gdzie

1

i

















2

1

n

i

in

1

n

i

k

k

n

h

e

e

2

e

u





















ale













2

sin

e

4

i

2

e

e

e

4

e

2

e

e

e

e

2

e

2

in

2

2

i

2

i

in

i

i

in

1

n

i

in

1

n

i











































































Po podstawieniu do









k

n

1

k

n

k

n

1

k

n

u

1

u

u

u





















mamy





2

sin

e

1

h

4

2

sin

e

h

4

e

e

2

in

k

2

2

in

1

k

2

in

k

in

1

k











































2

sin

h

4

1

2

sin

h

1

4

1

2

2

2

2























Po podzieleniu przez





in

k

e

otrzymujemy równanie:

Warunek konieczny stabilności Neumanna stwierdza, że schemat
różnicowy jest stabilny, jeżeli

1





Dla =0 mamy

2

sin

h

4

1

2

2











i na podstawie kryterium Neumanna otrzymujemy:

1

2

sin

h

4

1

2

2









czyli

0

2

sin

h

4

2

2

2













Dla = znajdujemy warunek na stosunek:

2

1

h

2





2

1

h

2





Warunkiem zbieżności schematu jawnego jest spełnienie

powyższego warunku. Warunek jest również prawdziwy

dla schematu jawnego w przypadku wielowymiarowego

równania parabolicznego.

Dowolne >0.

Ocenę prowadzimy przy =.





1

h

4

1

h

1

4

1

2

2

















czyli





2

2

2

h

4

1

h

1

4

1

h

4

1



























Prawa nierówność jest spełniona dla dowolnych , a z lewej

mamy:









4

h

2

1

2









4

h

2

1

2

Dla  spełniających nierówność:





1

h

4

1

h

1

4

1

2

2

















warunek:

jest spełniony dla dowolnego stosunku

2

h



W szczególności schemat Cranka-Nicholsona =0.5 jest stabilny

dla dowolnego stosunku kroków

2

h



Dla schematu o podwyższonej dokładności















12

h

5

.

0

2

mamy











4

h

2

1

2

bo





 4

h

12

h

2

2

Przedstawione rozważania można rozszerzyć na przypadki

wielowymiarowe jak również na równania o zmiennych

współczynnikach.

W przypadku równań wielowymiarowych ocena zbieżności

zależy również od sposobu aproksymacji warunków

brzegowych podobnie jak w przypadku równań eliptycznych.

 

 

 

 

 

 

 

 

t

e

L

,

t

u

t

e

0

,

t

u

x

w

x

,

0

,

u

x

w

x

,

0

u

,

u

v

,

u

L

P

1

t

0

xx

2

tt











Równania hiperboliczne

Jako przykład zostanie rozpatrzone równanie linii długiej
bez strat o długości L. Dla napięcia u mamy:

Wprowadzamy siatkę prostokątną:

x

t

0

N

h



i-1 i i+1

k+1
k
k-1

i funkcję węzłową oznaczamy:





ih

x

,

k

t

u

u

i

k

k

i









Przyjmujemy aproksymację pochodnych:

 

2

1

k

i

k

i

1

k

i

k

i

tt

v

u

u

2

u

Du

,

u















i

2

k

1

i

k

i

k

1

i

k

i

xx

h

u

u

2

u

u

,

u















i rozpatrujemy następujący schemat trójwarstwowy

z parametrem >0:









i

1

0

i

1

i

i

0

0

i

k

L

k

N

k

P

k

0

1

k

i

k

i

1

k

i

k

i

w

~

u

u

w

u

e

u

e

u

u

u

2

1

u

Du

































gdzie warunek początkowy

i

1

w

~ konstruujemy tak, aby

zachować rząd aproksymacji O(

2

).

Mamy

 

 

 

 

 

3

0

t

tt

2

0

t

t

O

t

,

x

,

u

5

.

0

t

,

x

,

u

0

,

x

u

,

x

u





















czyli

 

2

0

t

tt

i

0

t

t

i

0

i

1

i

O

,

u

5

.

0

,

u

u

u



















Z równania falowego mamy:

 

2

i

2

tt

i

h

O

u

v

,

u







Warunek początkowy dla pierwszej pochodnej będzie określony

z dokładnością O(h

2

+

2

), jeżeli przyjąć, że





2

2

i

0

2

i

1

0

i

1

i

h

O

w

v

5

.

0

w

u

u



















czyli

i

0

2

i

1

i

1

w

v

5

.

0

w

w

~









Ostatecznie schemat różnicowy dla rozwiązania równania

falowego jest









 

i

0

2

i

1

i

0

1

i

i

0

0

i

k

L

k

N

k

P

k

0

1

k

i

k

i

1

k

i

k

i

w

v

5

.

0

w

w

u

w

u

e

u

e

u

u

u

2

1

u

Du







































Ocena dokładności aproksymacji

Postępujemy podobnie jak poprzednio, a więc niech





k

i

k

i

k

i

t

,

x

u

u

z





gdzie

k

i

u

jest rozwiązaniem różnicowego zagadnienia









 

i

0

2

i

1

i

0

1

i

i

0

0

i

k

L

k

N

k

P

k

0

1

k

i

k

i

1

k

i

k

i

w

v

5

.

0

w

w

u

w

u

e

u

e

u

u

u

2

1

u

Du







































a u(x

i

,t

k

) jest rozwiązaniem problemu brzegowego:

 

 

 

 

 

 

 

 

t

e

L

,

t

u

t

e

0

,

t

u

x

w

x

,

0

,

u

x

w

x

,

0

u

,

u

v

,

u

L

P

1

t

0

xx

2

tt











w punkcie x

i

,t

k

.

Pisząc





k

i

k

i

k

i

t

,

x

u

z

u





otrzymujemy:









i

1

i

0

i

k

N

k

0

k

i

1

k

i

k

i

1

k

i

k

i

z

0

z

0

z

0

z

z

z

2

1

z

Dz



































gdzie



 

 



















 

i

0

2

i

1

i

i

0

i

k

i

1

k

i

k

i

1

k

i

k

i

w

v

5

.

0

w

1

,

x

u

w

t

,

x

Du

t

,

x

u

t

,

x

u

2

1

t

,

x

u









































Z konstrukcji warunku początkowego dla pochodnej wynika, że





2

2

i

h

O









Rozwijając w szereg Taylora mamy:



 







2

t

,

x

tt

k

i

1

k

i

1

k

i

k

i

,

u

t

,

x

u

2

t

,

x

u

t

,

x

u













Korzystając z otrzymanego wyniku mamy:



 









 







k

i

t

,

x

tt

2

k

i

k

i

1

k

i

1

k

i

,

u

t

,

x

u

t

,

x

u

2

1

t

,

x

u

t

,

x

u





















Stąd otrzymujemy z dokładnością do małych 4-go rzędu:



 

 











k

i

2

2

4

t

,

x

t

x

2

x

2

xx

1

k

i

k

i

1

k

i

,

u

,

u

12

h

,

u

t

,

x

u

t

,

x

u

2

1

t

,

x

u





























Podobnie





k

i

4

t

,

x

t

2

2

tt

2

k

i

,

u

v

12

,

u

v

1

t

,

x

Du







z dokładnością do małych 4-go rzędu. Ostatecznie otrzymujemy
ocenę błędu:

k

i

4

2

2

4

t

,

x

t

2

2

t

x

2

x

2

tt

2

xx

k

i

,

u

v

12

,

u

,

u

12

h

,

u

v

1

,

u

















Funkcja u spełnia równanie falowe, a więc

tt

2

xx

,

u

v

1

,

u 

stąd





2

2

k

i

h

O









niezależnie od 

!!!

Oznacza to również zależność odwrotną, a mianowicie

dobór  zależy tylko i wyłącznie od stabilności, a nie ma

wpływu na dokładność obliczeń.

Stabilność









 

n

0

2

n

1

n

0

1

n

n

0

0

n

k

N

k

0

1

k

n

k

n

1

k

n

k

n

w

v

5

.

0

w

w

u

w

u

0

u

0

u

u

u

2

1

u

Du







































Analizujemy stabilność schematu różnicowego przy jednorodnych
warunkach brzegowych:

Przyjmujemy rozwiązanie w postaci:











in

exp

u

k

k

n









1

k

n

k

n

1

k

n

k

n

u

u

2

1

u

Du





















i podstawiając do równania różnicowego:

po wykonaniu kilku przekształceń otrzymujemy:









2

sin

2

1

h

v

2

2

2

1

k

k

1

k

2

1

k

k

1

k





















































Dzieląc równanie przez 

k-1

i grupując wyrazy otrzymujemy

równanie określające :





0

1

2

sin

h

v

2

1

2

sin

h

v

2

2

1

2

1

1

2

2

2

2





















































0

1

r

4

1

r

2

1

2

2

2

2





























Badamy pierwiastki równania:

gdzie

h

v

r





przy =.

Wyróżnik:































2

2

2

2

r

1

4

1

r

4

1

r

4

Jeżeli <0, to równanie ma dwa sprzężone pierwiastki

zespolone 

1

, 

2

o module równym 1 co wynika ze wzoru

Viety:

1

2

1

2

1















Dla =0 otrzymujemy warunek Couranta:

1

h

v

r







który oznacza, że prędkość wędrówki fali na siatce h/

jest większa od prędkości fazowej.

1

h

v

r







Przypadek 0 prowadzi do pierwiastków większych co

do modułu od jedności i dlatego należy te przypadki

odrzucić.

Analizę można rozszerzyć na przypadki wielowymiarowe.

Wady i zalety metody różnicowej

Zalety:
1. Proste konstruowanie siatki podziałowej.
2. Prosta konstrukcja układu równań różnicowych szczególnie
w środowiskach izotropowych.
3. Opracowane oceny błędów metody i warunki stabilności.

Wady:
1. Duże trudności z dobrą aproksymacją brzegu lub wymuszony
mały krok siatki.
2. Trudności z utrzymaniem rzędu aproksymacji przy interpolacji
warunków brzegowych.
3. Praktycznie konieczność obliczania całego obszaru z tym
samym krokiem podziałowym.