Przykład

Równanie wahadła:

0

sin

2











Niech =1s

-2

Warunki początkowe:





1

.

2

0

t









około 86°





0

0

t







Sprowadzamy do układu równań I-go rzędu















sin





Warunki początkowe:

0

1

.

2

0

0











Obliczenia chcemy prowadzić z dokładnością 0.001

Startujemy z krokiem h=0.1. Krok wybrano jako

0.1 okresu wahadła liniowego.

Metoda Mersona

09972

.

0

3

K

sin

h

K

003324

.

0

3

K

h

K

09972

.

0

sin

h

K

0

h

K

1

0

2

1

0

2

0

1

0

1











































































099716

.

0

6

K

K

sin

h

K

003324

.

0

6

K

K

h

K

2

1

0

3

2

1

0

3































































099224

.

0

8

K

3

K

sin

h

K

0037394

.

0

8

K

3

K

h

K

3

1

0

4

3

1

0

4































































099701

.

2

K

4

K

3

K

sin

h

K

0098734

.

0

2

K

4

K

3

K

h

K

4

3

1

0

5

4

3

1

0

5







































































Błąd:

00013043

.

0

30

K

K

8

K

9

K

2

00032914

.

0

30

K

K

8

K

9

K

2

5

4

3

1

1

5

4

3

1

1





















































Dokładność założona została osiągnięta.
W następnym kroku można zwiększyć krok.

Rozwiązanie w chwili t=0.1

099386

.

6

K

K

4

K

49186

.

1

6

K

K

4

K

5

4

1

0

1

5

4

1

0

1











































i do następnego kroku możemy wystartować z nową
wartością kroku h

Metody włożone

lub

Metody Fehlberga – Runge -Kutty

Stosujemy metodę Runge – Kutty rzędu p i rzędu p+1
i aby zmniejszyć liczbę obliczanych współczynników
wybieramy je tak, że w obu metodach jest pierwszych
p współczynników K jednakowe, czyli





n

n

1

t

,

x

hf

K 

























h

c

t

,

K

a

x

hf

K

i

n

1

i

1

j

j

ij

n

i

i=2,3,..,p+1

i mamy dla metody rzędu p-go











p

1

i

i

p

i

n

1

n

K

w

x

x

a dla metody rzędu (p+1)-go















1

p

1

i

i

1

p

i

n

1

n

K

w

x

x

Ocenę błędu można zrobić stosunkowo prosto





 

2

p

1

p

n

n

p

p

1

i

i

p

i

n

1

n

h

O

h

t

,

x

K

w

x

x

























 

3

p

2

p

n

n

1

p

1

p

1

i

i

1

p

i

n

1

n

h

O

h

t

,

x

K

w

x

x



























Po odjęciu stronami otrzymujemy:









 

2

p

1

p

1

i

i

p

i

1

p

i

n

n

p

h

O

K

w

w

t

,

x



















gdzie

0

w

p

1

p





Znając błąd możemy postępować jak w metodzie
Mersona i rozwiązanie przyjmować z dokładniejszej
metody rzędu p+1.

Najczęściej stosowana metoda RKF45 ma współczynniki



























































h

13

12

t

,

2197

K

7296

K

7200

K

1932

x

hf

K

h

8

3

t

,

32

K

9

K

3

x

hf

K

h

25

.

0

t

,

K

25

.

0

x

hf

K

t

,

x

hf

K

n

3

2

1

n

4

n

2

1

n

3

n

1

n

2

n

n

1



















































2

h

t

,

K

40

11

K

4104

1859

K

2565

3544

K

2

K

27

8

x

hf

K

h

t

,

K

4104

845

K

513

3680

K

8

K

216

439

x

hf

K

n

5

4

3

2

1

n

6

n

4

3

2

1

n

5

6

5

4

3

1

n

K

55

2

K

50

9

5

1

K

4104

2197

56430

28561

K

2565

1408

12825

6656

K

216

25

135

16













 























































Błąd

Rozwiązanie wykorzystując metodę dokładniejszą jest

6

5

3

2

1

n

1

n

K

55

2

K

50

9

K

56430

28561

K

12825

6656

K

135

16

x

x

















Metoda gwarantuje obliczenia z błędem rzędu h

4

.

Stabilność metod Rungego - Kutty

Przedziały stabilności absolutnej dla metody Rungego-Kutty
wynoszą:

Ax

x 





Dla metody rzędu pierwszego:

2

h

0







rzędu drugiego:

2

h

0







rzędu trzeciego:

51

.

2

h

0







rzędu czwartego:

78

.

2

h

0







dla każdej wartości własnej macierzy A.

Żadna z metod Rungego-Kutty czy Fehlberga-Rungego-Kutty

nie jest odpowiednia do równań sztywnych

Typowe równania różniczkowe cząstkowe rzędu drugiego

Przykłady fizyczne

Do przewodnika o przewodności przyłożono napięcie U.

Znaleźć rozkład gęstości prądu w przewodniku.

x

y

z

j

y

(x,y,z)

j

y

(x,y+y,z)

y

x

z

j

x

(x,y,z)

j

x

(x+x,y,z)

j

z

(x,y,z)

j

z

(x,y,z+z)

x

y

z

j

y

(x,y,z)

j

y

(x,y+y,z)

y

x

z

j

x

(x,y,z)

j

x

(x+x,y,z)

j

z

(x,y,z)

j

z

(x,y,z+z)

Bilans:

























0

y

x

z

z

,

y

,

x

j

y

x

z

,

y

,

x

j

z

x

z

,

y

y

,

x

j

z

x

z

,

y

,

x

j

z

y

z

,

y

,

x

x

j

z

y

z

,

y

,

x

j

y

y

y

y

x

x













































































0

y

x

z

z

,

y

,

x

j

y

x

z

,

y

,

x

j

z

x

z

,

y

y

,

x

j

z

x

z

,

y

,

x

j

z

y

z

,

y

,

x

x

j

z

y

z

,

y

,

x

j

y

y

y

y

x

x





















































Rozwijamy w szereg Taylora w otoczeniu punktu (x,y,z):





 

 

2

z

,

y

,

x

x

x

x

x

O

x

x

j

)

z

,

y

,

x

(

j

)

z

,

y

,

x

x

(

j























 

 

2

z

,

y

,

x

y

y

y

y

O

y

y

j

)

z

,

y

,

x

(

j

)

z

,

y

y

,

x

(

j























 

 

2

z

,

y

,

x

z

z

z

z

O

z

z

j

)

z

,

y

,

x

(

j

)

z

z

,

y

,

x

(

j



















Podstawiając otrzymujemy:

0

V

z

j

y

j

x

j

z

y

x

































gdzie

z

y

x

V











ostatecznie:

0

z

j

y

j

x

j

z

y

x



















Wektor gęstości prądu j wyrażamy przez potencjał :











j

czyli

z

j

y

j

x

j

z

y

x





































lub

0





 j

Podstawiając

do równania:

0

z

j

y

j

x

j

z

y

x



















z

j

y

j

x

j

z

y

x





































mamy:

0

z

z

y

y

x

x















































































Jeżeli przewodność ośrodka  jest stała, to

0

y

y

x

2

2

2

2

2

2

























Równanie nazywamy równaniem Laplace’a (eliptycznym)

Prosty przykład:

Dana jest przewodząca płytka prostokątna o wymiarach 2ax2b
i stałej grubości h. Przewodność elektryczna płytki wynosi 

Na jednym z boków z o długości 2a znajduje się elektroda,
której potencjał jest pokazany na rysunku. Druga elektroda
znajdująca się na przeciwległym boku jest uziemiona.
Wyznaczyć rozkład gęstości prądu w płytce i moc traconą w niej.

2b

2a

h

-a

a

l

U(l)

E

1. FIZYKA

Jak sobie wyobrażamy rozpływ gęstości prądu w płytce?

Jakie stawiamy założenia upraszczające?

Jakie są warunki zadania?

Na podstawie tych rozważań budujemy

model matematyczny

W naszym przypadku mamy następujące wnioski:
1. Można przyjąć dwuwymiarowy model matematyczny

2D

2. Układ współrzędnych prostokątnych (x,y)

Jak umieścić układ
współrzędnych?

a

l

U(l)

U

m

-a

x

y

-a

a

0

2b

Model matematyczny

Wektor gęstości prądu
j ma dwie składowe j

x

, j

y

będące funkcjami x i y.
Spełnia równanie:

0





 j

jest związany z natężeniem
pola elektrycznego E:

E

j 



a pole elektryczne spełnia równania:

0

0













E

E

ponieważ materiał jest jednorodny i izotropowy, więc
równania:

0





 E

i

 

0























E

E

j

są

równoważne

. Wystarczy określić rozkład pola elektrycznego

a z prawa Ohma wyznaczymy rozkład gęstości prądu.
Zadanie sprowadza się więc do rozwiązania układu równań:

0

0













E

E

Ze względu na pierwsze przyjmujemy:







E

i podstawiając do drugiego równania mamy:

0

2







lub

0

y

x

2

2

2

2

















x

-a

a

0

2b

a

l

U(l)

U

m

Symetrie i warunki brzegowe:

j

x

(x=0,y)=?

-a

x

-a

a

0

2b

a

l

U(l)

U

m

-a

0

j

0

x

x





ale

x

E

j

x

x

















czyli

0

x

0

x











i co więcej









y

,

x

y

,

x









x

-a

a

0

2b

0

j

0

x

x





Fizyka:

a

x

0

x











a

x

x

0

j







0

y

0







uziemiona elektroda

b

2

y

m

a

x

a

U









Ostatecznie model matematyczny ma postać:

0

y

x

2

2

2

2

















a

x

0

x











0

x

0

x











0

y

0







b

2

y

m

a

x

a

U









Przedstawiamy potencjał w postaci:





   

y

Y

x

X

y

,

x





i podstawiając mamy:

0

dy

Y

d

X

Y

dx

X

d

2

2

2

2





Dzieląc przez XY mamy:

0

dy

Y

d

Y

1

dx

X

d

X

1

2

2

2

2





czyli:

2

2

2

dy

Y

d

Y

1





i

2

2

2

dx

X

d

X

1







Mamy:

0

Y

dy

Y

d

0

X

dx

X

d

2

2

2

2

2

2













i rozwiązanie pierwszego równania jest:

 

 

 

x

sin

C

x

cos

C

x

X

2

1

















y

,

x

y

,

x









ale

stąd

0

C

2



czyli

 

 

x

cos

C

x

X

1





Z drugiego warunku:

a

x

0

x











mamy:

 

0

a

sin

C

1





ponieważ

0

C

1



więc

 

0

a

sin





czyli

a

n

n







czyli

 

















 



0

n

n

1

a

x

n

cos

C

x

X

0

Y

dy

Y

d

2

n

2

2







Drugie równanie:

ma rozwiązanie:

 









y

sinh

C

y

cosh

C

y

Y

n

n

4

n

n

3

n









Z warunku brzegowego:

0

y

0







mamy:

0

C

n

3



Biorąc pod uwagę:





   

y

Y

x

X

y

,

x





mamy:



























1

n

n

n

n

y

sinh

x

cos

C

y

,

x

i pisząc C

n

=C

1n

C

3n

Z ostatniego warunku brzegowego:

b

2

y

m

a

x

a

U









mamy:























1

n

m

n

n

n

a

x

a

U

b

2

sinh

x

cos

C

Korzystając z ortogonalności funkcji cos liczymy współczynniki C

n















 















 











 

a

0

m

a

0

2

n

dx

x

a

n

cos

a

x

a

U

a

b

n

2

sinh

dx

x

a

n

cos

C

po wykonaniu całkowania
mamy:







































a

b

1

n

2

2

sinh

1

1

n

2

2

C

0

C

2

1

n

2

n

2















































































a

x

1

n

2

cos

a

b

1

n

2

2

sinh

a

y

1

n

2

sinh

1

n

2

2

U

y

,

x

0

n

2

m

Znając potencjał możemy określić rozkład gęstości prądu
z równania:











j

czyli:

y

j

x

j

y

x































































































a

x

1

n

2

sin

a

b

1

n

2

2

sinh

a

y

1

n

2

sinh

1

n

2

4

a

U

y

,

x

j

0

n

m

x









































































a

x

1

n

2

cos

a

b

1

n

2

2

sinh

a

y

1

n

2

cosh

1

n

2

4

a

U

y

,

x

j

0

n

m

x

Obliczenie mocy traconej w płytce







V

2

dv

P

j

Uwzględniając warunki zadania mamy:













b

2

0

a

0

2

y

2

x

dy

dx

j

j

h

2

P

Podstawiając i wykonując całkowanie otrzymujemy:

























































0

n

2

2

m

a

b

1

n

2

2

sinh

1

n

2

4

U

a

hb

2

P

Przykład

-h

h

d

-d

x

y

Do prostokątnej płytki o wymiarach 2hx2d
i stałej grubości H przyłożono dwie elektrody.

2g

Górna elektroda położna w środku

o szerokości 2g i potencjale V i dolna

elektroda wzdłuż dolnego boku

o potencjale 0.

Przewodność płytki jest stała

i wynosi .

Wyznaczyć rozkład gęstości prądu
w płytce i rezystancję zastępczą płytki
przy tak przyłączonych elektrodach.

Opis matematyczny:

E

j

j

E

















0

0

Przyjmując:







E

mamy:











j

a potencjał 

spełnia równanie Laplace’a:

0

2







-h

h

d

-d

x

y

2g

Wniosek z geometrii elektrod:

Potencjał  jest jedynie funkcją

współrzędnych x,y.

Potencjał  jest funkcją parzystą zmiennej x

czyli









y

,

x

y

,

x









co oznacza, że jako rozwiązanie należy przyjąć funkcję
parzysta względem x i można nasze zadanie rozważyć
w obszarze

d

y

d

h

x

0











-h

h

d

-d

x

y

2g

Warunki brzegowe:

d

y

h

x

0















0

d

,

x







d

y

d

0

x











0

x

0

x











i

d

y

d

h

x











h

x

0

x











-h

h

d

-d

x

y

2g

i ostatni:

  



 





 

d

y

V

x

y

x

1

y

,

x

x























gdzie

 



















h

x

g

dla

0

g

x

0

dla

1

x





   

y

Y

x

X

y

,

x





Jak poprzednio przyjmujemy:

i mamy:

0

dy

Y

d

Y

1

dx

X

d

X

1

2

2

2

2





2

2

2

dy

Y

d

Y

1





2

2

2

dx

X

d

X

1







oraz

które mają rozwiązanie:

 

 

 

x

sin

C

x

cos

C

x

X

2

1









 

 

 

y

sinh

C

y

cosh

C

y

Y

4

3









Z warunku:

0

x

0

x











wynika

0

C

2





czyli

0

C

2



i stąd

 

 

x

cos

C

x

X

1





h

x

0

x











Z warunku brzegowego:

wynika:

 

0

h

sin

C

1





Ponieważ

0

C

1



więc

h

n

n







Z warunku brzegowego:





0

d

,

x







mamy:

 

 

 

0

d

sinh

C

d

cosh

C

d

Y

4

3













czyli

 

d

ctgh

C

C

3

4





Potencjał będzie po wprowadzeniu zastępczej stałej:





















x

cos

d

sinh

d

y

sinh

C

y

,

x

n

1

n

n

n

n



















gdzie

h

n

n







Ostatni warunek brzegowy:

  



 





 

d

y

V

x

y

x

1

y

,

x

x























daje:

 





 













 

















































1

n

n

n

n

n

n

n

V

x

x

cos

d

sinh

d

2

cosh

x

1

d

cosh

x

2

C

no i mamy schody!!!

Jak wybrnąć

z tych kłopotów?

Dla skrócenia zapisu wprowadzamy oznaczenie:

 





 









 

































h

,

g

x

dla

d

sh

d

2

ch

g

,

0

x

dla

d

ch

2

x

n

n

n

1

n

n

0

n

n

czyli

  



 















1

n

n

n

n

V

x

x

cos

x

C

Żądamy minimum błędu aproksymacji średniokwadratowej:





  



 

min

dx

V

x

x

cos

x

C

,...

C

,...,

C

,

C

F

h

0

2

1

n

n

n

n

n

2

1

























 





Pozostała już tylko arytmetyka!

Obliczamy ekstremum funkcji wielu zmiennych:

0

dC

dF

k



i mamy:

  



 

  



0

dx

x

cos

x

V

x

x

cos

x

C

h

0

k

k

1

n

n

n

n



























 





k=1,2,...

Otrzymujemy nieskończony układ liniowych równań:



 





 

























































g

0

k

1

n

h

g

k

n

1

k

1

n

g

0

k

n

0

k

0

n

n

dx

x

cos

V

dx

x

cos

x

cos

dx

x

cos

x

cos

C

Niestety całki:



 





 



0

dx

x

cos

x

cos

0

dx

x

cos

x

cos

h

g

k

n

g

0

k

n

















i układ równań ma nieskończoną liczbę niewiadomych.

Rozwiązujemy w ten sposób, że ograniczamy liczbę wyrazów

i rozwiązujemy układ o skończonej liczbie niewiadomych.

Jest to jednak metoda bardzo pracochłonna i wymagająca

albo dobrej znajomości metod rozwiązywania równań

o nieograniczonej liczbie niewiadomych

albo kilkukrotnego rozwiązania odpowiednio powiększanej

liczby równań i ocenie odrzuconej części.

W takiej sytuacji bezwzględnie bardziej efektywne są metody

numeryczne.

Przejdziemy do omówienia metod numerycznych stosowanych
do rozwiązywania zagadnień opisanych równaniem eliptycznym.
Jest to następująca grupa zagadnień:

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.