Transcript MN7-2012

METODY NUMERYCZNE
Wykład 7
Równania różniczkowe - przegląd
dr hab.inż. Katarzyna Zakrzewska, prof.AGH,
Katedra Elektroniki, AGH
e-mail: [email protected]
http://home.agh.edu.pl/~zak
1
Równania różniczkowe - wprowadzenie
Równania różniczkowe są popularnie spotykane we
wszystkich dziedzinach nauk ścisłych i przyrodniczych a
szczególnie w:
• Fizyce (np. równania Maxwell’a)
• Mechanice (np. równania ruchu harmonicznego)
• Elektronice (np. stany nieustalone w obwodach
elektrycznych)
• Automatyce (np. warunki sterowalności układu)
• i wielu innych dziedzinach nauki i techniki
2
Równania różniczkowe zwyczajne
Równania różniczkowe zwyczajne – jest to równanie w którym
występują stałe oraz funkcje niewiadome i pochodne funkcji
niewiadomych zależne od jednej zmiennej niezależnej.
n
d y
dx
n
 kn
d
n 1
dx
y
n 1
2
 .........  k 3
Przykład:
M
dt
2
dx
2
 k2
dy
dx
b
dx
dt
 k1 y  F ( x )
b
K
2
d x
d y
 Kx  0
M
x
3
Cząstkowe równania różniczkowe
Cząstkowe równanie różniczkowe – jest to równanie zawierające
funkcję niewiadomą dwóch lub więcej zmiennych oraz
niektóre z jej pochodnych cząstkowych.
ut  b   u  0
Jednym z najprostszych równań
różniczkowych cząstkowych
jest równanie transportu:
u  u (t , x )
 u
u
 u  
,...,
xn
  x1




4
Cząstkowe równania różniczkowe
Zauważamy, że pochodna kierunkowa funkcji u w kierunku
wektora v=(1,b) є Rn+1 znika. Zatem ustalając dowolny
punkt (t,x) є R+ x Rn i kładąc dla s є R dostajemy:
z ( s )  u ( t  s , x  sb )
dz ( s )
ds
 u t ( t  s , x  sb )   x u ( t  s , x  sb )  0
Zatem z(s) jest funkcją stałą.
Ustalając wartość rozwiązania na każdej prostej równoległej do
wektora (1,b) dostajemy rozwiązanie zadania.
5
Cząstkowe równania różniczkowe – zagadnienia
początkowe
Zagadnienia początkowe, zakładamy, że w chwili t=0 zadana
jest wartość funkcji u(0,x). Wówczas zagadnienie
początkowe:
 u  b   u  0 dla t  0 , x  R
 t

n
 u ( 0 , x )  g ( x ) dla x  R
n
ma rozwiązanie:
u ( t , x )  g ( x  tb ) dla t  0 , x  R
n
Jeśli funkcja g jest klasy C1 to rozwiązanie równania jest
rozwiązaniem klasycznym oraz jest ono jednoznaczne
6
Rozwiązywanie zagadnień
początkowych
Wprowadzimy teraz kilka oznaczeń, niech:
 Y(x) – oznacza dokładne rozwiązanie
 y(x) – oznacza rozwiązanie przybliżone
dY ( x i )
Y i  Y ( x i ),
Yi 
y i  y ( x i ),
y  f ( xi , yi )  f i
'
dx
 f ( x i , Yi )
'
i
gdzie :
x i  ( x 0 ; b , i  1, 2 ,...., N
są punktami, w których
wyznaczamy przybliżone
rozwiązania
7
Błąd metody
Wielkość Tn nazywamy błędem metody powstałym
przy przejściu od xn do xn+1
x i  x 0  ih
h – krok całkowania
Błąd metody możemy wyrazić jako funkcję zmiennej h
i przedstawić w postaci:
Tn  h
p 1

Gdzie stała
p
 O (h

p
p2
)
 0 , to liczbę p będziemy nazywać
rzędem metody przybliżonej
8
Cząstkowe równania różniczkowe – zagadnienia
niejednorodne
W celu rozwiązania zagadnienia niejednorodnego:
 u  b   u  f dla t  0 , x  R
 t

n
 u ( 0 , x )  g ( x ) dla x  R
podstawmy:
wówczas:
n
z ( s )  u ( t  s , x  sb )
dz ( s )
ds
 u t ( t  s , x  bs )   u ( t  s , x  bs )  b 
 f ( t  s , x  bs )
9
Cząstkowe równania różniczkowe – zagadnienia
niejednorodne
Zatem:
y ( t , x )  g ( x  tb )  z ( 0 )  z (  t ) 
0
dz
t
ds


0


f ( t  s , x  sb ) ds 
t
t


f ( s , x  ( s  t ) b ) ds
0
Rozwiązaniem zagadnienia jest więc:
u ( t , x )  g ( x  tb ) 

t
f ( s , x  ( s  t ) b ) ds
0
10
Całka zupełna dla równań rzędu 1
Ogólne równanie różniczkowe cząstkowe rzędu 1 można
zapisać w postaci:
F ( x, u ,  u )  0
gdzie:
x  , u :   R
 u    x 1 u ,...,  x 2 u 
F : R R  R
n
11
Całka zupełna dla równań rzędu 1 – zagadnienia
brzegowe
Ogólne równanie różniczkowe cząstkowe rzędu 1 można
zapisać w postaci, spróbujemy rozwiązać warunki brzegowe
(zagadnienie Cauchy’ego)
F ( x , u ,  u )  0 w obszarze
R
n
u  g na powierzchn i    
Zakładamy, że funkcje F i g oraz powierzchnia Г są dostatecznie gładkie.
12
Całka zupełna dla równań rzędu 1 – zagadnienia
brzegowe
Staramy się połączyć punkt xєΩ z pewnym punktem x0єГ pewną
krzywą ɣ w taki sposób, aby można było policzyć wartości
rozwiązania u wzdłuż tej krzywej:
  x ( s ), s  I  R 
z ( s )  u ( x ( s )), p ( s )   u ( x ( s ))
Chcemy dobrać krzywą x(s) tak, aby można było policzyć z(s) i
p(s). W tym celu policzymy dp(s)/d(s):
p (s) 
i
u
j 1
 u
2
n
( x ( s ))  x ( s ) dla i  1,..., n u ij 
j
ij
xi x j
13
Całka zupełna dla równań rzędu 1 – zagadnienia
brzegowe
Z drugiej strony różniczkując równanie ogólne względem xi:
F
xi
( x, u ,  u ) 
F
z
n
( x , u ,  u )u xi 
F
 p
j 1
( x , u ,  u ) u ij  0 i  1,... n
j
Zakładamy, że:
x (s) 
j
F
p j
( x ( s ), z ( s ), p ( s )) dla j  1,..., n
14
Całka zupełna dla równań rzędu 1 – zagadnienia
brzegowe
Otrzymujemy wtedy:
p (s) 
i
F
xi
( x ( s ), z ( s ), p ( s )) 
F
z
 x ( s ), z ( s ),
p ( s )  p ( s ) dla i  1,..., n
i
Ostatecznie otrzymujemy:

z(s) 
n
u
 x
j 1
j
n
( x ( s )) x ( s ) 
j

j 1
j
p (s)
F
p j
( x ( s ), z ( s ), p ( s ))
15
Całka zupełna dla równań rzędu 1 – zagadnienia
brzegowe
Stosując zapis wektorowy otrzymujemy układ (2n+1) równań
różniczkowych zwyczajnych zwany układem charakterystyk
całki zupełnej.
p ( s )   x F ( x ( s ), z ( s ), p ( s ))   z F ( x ( s ), z ( s ), p ( s ))  p ( s )
z ( s )   p F ( x ( s ), z ( s ), p ( s ))  p ( s )
x ( s )   p F ( x ( s ), z ( s ), p ( s ))
16
Całka zupełna dla równań rzędu 1 – zagadnienia
brzegowe - przykład
Rozpatrzmy układ:
  x 1  0 , x 2  0 
 x1 u x 2  x 2 u x 1  u

 u ( x1 , x 2 )  g ( x1 )
  x1  0 , x 2  0 
Wówczas równania charakterystyk mają postać:

x1   x 2 ,


x 2  x1 , z  z
17
Całka zupełna dla równań rzędu 1 – zagadnienia
brzegowe - przykład
Rozwiązując ten układ równań z uwzględnieniem warunku
brzegowego dostajemy
x1 ( s )  x 0 cos s , x 2 ( s )  x 0 sin s , z ( s )  g ( x 0 ) e
x0  0,
0 s
s

2
Ostatecznie rugując parametr s dostajemy rozwiązanie:

x2 
u ( x1 , x 2 )  g ( x  x ) exp  arctg
 , x1  1, x 2  0
x1 

2
1
2
2
18
Równanie liniowe rzędu 2
Ogólne równanie różniczkowe cząstkowe liniowe rzędu 2
określone w obszarze Ω ↄ Rn ma postać:
n
a
n
k ,l
( x ) u xk , xl ( x ) 
k , l 1
b
k
( x ) u xk ( x )  c ( x ) u  f ( x )
k 1
Ogólne równanie cząstkowe drugiego rzędu dwóch zmiennych
niezależnych liniowe możemy zapisać:
 f
2
A( x, y )
x
2
 f
2
 2 B ( x, y )
 xy
 f
2
 C ( x, y )
y
2
 F ( x, y, f ,
f
,
f
x y
A,B,C są określone w pewnym obszarze Ω ↄ R2 a F jest
określona na Ω ↄ R3
19
)
Transformacja Laplace’a
Całką Laplace’a funkcji f nazywamy całkę niewłaściwą
postaci:

e
 st
f ( t ) dt  F ( s )
0
Przy czym o funkcji f(t) zakładamy że jest funkcją
rzeczywistą zmiennej rzeczywistej t określoną dla
każdej wartości t > 0 i przedziałami ciągłą.
Ponieważ całka jest całką niewłaściwą to nie dla
wszystkich funkcji f(t) spełniających podane warunki
jest ona zbieżna.
20
Transformacja Laplace’a c.d.
Funkcję f(t) dla których istnieje całka Laplace’a
nazywamy oryginałami natomiast odpowiadające i
funkcje F(s) nazywamy transformatami Laplace’a.
Samą transformację Laplace’a zwaną także
przekształceniem Laplace’a w środowisku
inżynierskim często zapisujemy jako:

F ( s )  L { f ( t )} 
e
 st
f ( t ) dt
0
21
Transformacja odwrotna Laplace’a
Transformacja odwrotna Laplace’a dla klasy funkcji
spełniających podane założenia wyraża się
wzorem:
1
L { F ( s )}  f ( x ) 
1
2 i
x  iu
lim
u  

st
F ( s )e ds
x  iu
22
Transformacja Laplace’a – przykładowe
funkcje
x  0 F ( s )  L{ f ( x )}
f ( x ),
e

1
1
x
Transformata pochodnej:
L { f ( x )} 
n
s
s
n 1
 s F (s)  f
n
1
ax
 sx
f ( x ) dx 
n
0
n!
n
e
n 1
( 0 )  ........  s
n 1
f (0)
sa
a
sin ax
cos ax
sinh ax
cosh ax
s a
2
2
s
s a
2
t
2
a
s a
2
Transformata całki:
2
L {  f ( )d  } 
0
1
F (s)
s
s
s a
2
2
23
Własności transformacji Laplace’a
• Linowość:
L{ af ( x )  bg ( x )  ch ( x )}  aL ( f ( x ))  bL ( g ( x ))  cL ( h ( x ))
gdzie a, b, c to współczynniki
• Przesunięcie w dziedzinie zmiennej zespolonej:
Jeśli
L { f ( x )}  F ( s )
to
L{e f ( x )}  F ( s  a )
to
s
L { f ( ax )}  F  
a a
at
• Zmiana skali
Jeśli
L { f ( x )}  F ( s )
1
24
Transformacja Laplace’a - przykład
Rozwiązywanie równania różniczkowego przy pomocy
transformacji Laplace’a:
3
dy
 2y  e
x
y (0)  5
dx
Krok 1: Transformacja obydwu stron równania różniczkowego
 dy

x
L 3
 2 y   L e 
 dx

3[ sY ( s )  y ( 0 )]  2 Y ( s ) 
1
s 1
Uwzględniając warunek początkowy y(0) = 5 otrzymujemy:
3[ sY ( s )  5 ]  2 Y ( s ) 
1
s 1
25
Transformacja Laplace’a – przykład c.d.
Krok 2: Wyrażanie Y(s) jako funkcji s
( 3 s  2 )Y ( s ) 
( 3 s  2 )Y ( s ) 
Y (s) 
1
s 1
 15
15 s  16
s 1
15 s  16
( s  1)( 3 s  2 )
Zapisanie Y(s) w postaci ułamków prostych
15 s  16
( s  1)( 3 s  2 )
Y (s) 
1
s 1


A
s 1

B
3s  2
A   1;
B  18 ;
18
3s  2
26
Transformacja Laplace’a – przykład c.d.
Krok 3: Odwrotna transformacja obydwu stron równania
Y (s) 
1
s 1

18
3s  2

1
s 1

6
s  0 . 666667
6
 1 

1 
L {Y ( s )}  L 
 L 

 s 1
 s  0 . 666667 
1
Uwzględniając że:
1
 1 
 at
L 
e
sa
1
Rozwiązanie równania wynosi:
y( x)  e
x
 6e
 0 . 666667 x
27
Równanie Poissona
Równaniem Poissona nazywamy niejednorodne
równanie Laplace’a
u (r , t )  f (r , t )
r  x , y , z  R
 

2
x
2


3
2
y
2


2
z
2
28
Równanie Poissona
Równanie Poissona możemy podać explicite dla
przestrzeni 2 i 3 wymiarowej:

2
x

2
u ( x, y ) 

y
2
x
2
2
u ( x, y, z ) 
2

u ( x, y )  f ( x, y )
2
y
2
u ( x, y, z ) 

2
z
2
u ( x, y, z )  f ( x, y, z )
29
Równanie Poissona
Rozwiązanie równania Poissona wyrażamy za pomocą
funkcji Greena:
u (r ) 
 G ( r , r ' ) f ( r ' ) dr '
R
G (r , r ' ) 
1
4 r  r '
30
Rozwiązywanie równań różniczkowych
zwyczajnych – metoda Eulera
Metoda Eulera pozwala na numeryczne rozwiązanie
równania różniczkowego zwyczajnego rzędu
pierwszego postaci:
dy
dx
 f  x , y ,
y 0   y 0
Przykłady:
dy
 2 y  1 .3e
x
,
dx
dy
 1 .3e
x
 2 y,
dx
f  x , y   1 .3e
x
y 0   5
e
y 0   5
dy
 2y
y
dy
y 0   5
 x y  2 sin( 3 x ),
2
2
dx
2 sin( 3 x )  x y
2

dx
f x, y  
e
y
2
,
2 sin( 3 x )  x y
2
e
y 0   5
2
y
31
Rozwiązywanie równań różniczkowych
zwyczajnych – metoda Eulera
Dla x = 0 wartość y = y0 przyjmując że x = x0 = 0
y
wartość prawdziwa
( x0 , y0 )
y1, wartość
przewidywana
Φ
krok h
x1
x
Znając f(x, y) i mając dane wartości x0 i y0 z
warunku początkowego y(x0) = y0 można obliczyć
nachylenie funkcji f(x, y) do osi X w punkcie (x0, y0)
32
Rozwiązywanie równań różniczkowych
zwyczajnych – metoda Eulera
tan  
y1  y 0
x1  x 0
 f ( x0 , y0 )
Po przekształceniu otrzymujemy:
y 1  y 0  f  x 0 , y 0  x1  x 0 
Oznaczając x1-x0 jako krok h otrzymujemy:
y 1  y 0  f  x 0 , y 0 h
Wykorzystując obliczaną wartość y1 można obliczyć
wartość y2 dla x2 jako:
y 2  y 1  f  x 1 , y 1 h
x 2  x1  h
33
Rozwiązywanie równań różniczkowych
zwyczajnych – metoda Eulera
Można zatem wyprowadzić wzór rekurencyjny:
y i 1  y i  f  x i , y i h
y
wartość prawdziwa
yi+1, wartość przewidywana
Φ
yi
h
krok
xi
xi+1
x
34
Metoda Eulera - Przykład
Kula nagrzana do temperatury 1200 K schładza się do
temperatury 300K. Zakładając że w procesie
schładzania ciepło jest oddawane do otoczenia
jedynie przez radiację to temperatura kuli opisana
jest równaniem:
d
dt
  2 . 2067  10
 12

4

 81  10 ,
8
 0   1200 K
Jaka jest temperatura kuli po 480 sekundach jej
schładzania?
f t ,     2 . 2067  10
 12

4
 81  10
8

35
Metoda Eulera – Przykład c.d.
Wzór rekurencyjny metody Eulera:  i 1   i  f t i ,  i h
Zakładamy krok h = 240
Dla i = 0, t0 = 0, Θ0 = 1200:
 1   0  f t 0 ,  0 h  1200  f 0 ,1200   240

 1200   2 . 2067  10
 12
1200
4
 81  10
8
  240
 1200    4 . 5579   240  106 . 09 K
t  t1  t 0  h  0  240  240
36
Metoda Eulera – Przykład c.d.
Dla i = 1, t1 = 240, Θ0 = 106.09:
 2   1  f t1 ,  1 h  106 . 09  f  240 ,106 . 09   240

 106 . 09   2 . 2067  10
 106 . 09  0 . 017595
 12
  240
106 . 09
4
 81  10
8
  240
 110 . 32 K
t  t 2  t1  h  240  240  480
Po wykonaniu dwóch iteracji metody Eulera
otrzymujemy że temperatura kuli po 480 s
wyniesie 110.32K
Czy to prawda?
37
Metoda Eulera – Przykład c.d.
temperatura Θ(K)
Porównanie rozwiązania dokładnego z rozwiązaniem
otrzymanym przy użyciu metody Eulera.
dokładne rozwiązanie
czas t(s)
38
Metoda Eulera – Przykład c.d.
Dobór odpowiedniego kroku h jest kluczowy dla
odpowiedniej dokładności rozwiązania stosując
metodę Eulera
temperatura Θ(K)
dokładne rozwiązanie
czas t(s)
rozmiar
kroku h   480
480
240
120
60
30

| t | %
-987.81
110.32
546.77
614.97
632.77
252.54
82.964
15.566
5.0352
2.2864
39
Rozwiązywanie równań różniczkowych
zwyczajnych – metoda Rungego - Kutty
Metoda Rungego-Kutty pozwala na numeryczne
rozwiązanie równania różniczkowego zwyczajnego
pierwszego rzędu postaci:
dy
 f  x , y ,
y 0   y 0
dx
Korzystając z rozwinięcie w szereg Taylora:
y i 1  y i 
dy
dx
 x i 1
 xi  
xi , yi
 y i  f ( x i , y i )  x i 1  x i  
1
2!
2
1 d y
2! dx
 x i 1
2
 xi  
3
1 d y
2
xi , yi
f ' ( x i , y i )  x i 1  x i  
2
3! dx
1
3!
 x i 1
3
 x i   ...
3
xi , yi
f ' ' ( x i , y i )  x i  1  x i   ...
3
Widać że metoda Eulera jest metodą Rngego-Kutty pierwszego rzędu
y i 1  y i  f  x i , y i h
40
Rozwiązywanie równań różniczkowych
zwyczajnych – metoda Rungego - Kutty
Wzór dla metody Rungego-Kutty drugiego rzędu
będzie wyglądał następująco:
y i  1  y i  f  x i , y i h 
1
2!
2
f  x i , y i h
dla metody czwartego rzędu:
y i  1  y i  f  x i , y i h 
1
2!
f
'
 x i , y i h
2

1
3!
f
''
 x i , y i h
3

1
4!
f
'''
 x i , y i h 4
Jak wyznaczyć pochodne f’ metody stopnia drugiego i
f’, f’’, f’’’ dla metody stopnia czwartego?
41
Rozwiązywanie równań różniczkowych
zwyczajnych – metoda Rungego - Kutty
Dla metody Rungego-Kutty drugiego rzędu wzór:
y i  1  y i  f  x i , y i h 
można zapisać jako:
gdzie:
k1  f xi , y i 
1
2!
2
f  x i , y i h
y i 1  y i  a 1 k 1  a 2 k 2 h
k 2  f  x i  p 1 h , y i  q 11 k 1 h 
aby wyznaczyć współczynniki a1, a2, p1, q11 należy
rozwiązać kład równań:
1
1
a 2 p1 
a 2 q 11 
a1  a 2  1
2
2
zazwyczaj wartość a2 wybiera się aby rozwiązać pozostałe
42
Rozwiązywanie równań różniczkowych
zwyczajnych – metoda Rungego - Kutty
Zazwyczaj a2 przyjmuje jedną z trzech wartości: ½, 1, 2/3
Metoda Heun’a
a2 
a1 
1
2
Metoda midpoint
1
a2  1
2
p1  1
Metoda Raltson’a
q 11  1
a1  0
p1 
1
2
q 11 
1
2
a1 
1
3
a2 
2
p1 
3
3
4
q 11 
3
4
1
1

y i 1  y i   k 1  k 2  h
2
2

y i 1  y i  k 2 h
2
1

y i 1  y i   k 1  k 2  h
3
3

k1  f xi , y i 
k1  f xi , y i 
k1  f xi , y i 
k 2  f  x i  h, y i  k1 h 
1
1
3
3




k 2  f  xi  h, y i  k1h  k 2  f  xi  h, y i  k1h 
2
2
4
4




43
Metoda Rungego - Kutty - Przykład
Kula nagrzana do temperatury 1200 K i schładza się do
temperatury 300K. Zakładając że w procesie
schładzania ciepło jest oddawane do otoczenia
jedynie przez radiację to temperatura kuli opisana
jest równaniem:
d
dt
  2 . 2067  10
 12

4

 81  10 ,
8
 0   1200 K
Jaka jest temperatura kuli po 480 sekundach jej
schładzania?
f t ,     2 . 2067  10
 12

4
 81  10
8

44
Metoda Rungego - Kutty – Przykład c.d.
Dla metody Heun’a:
1
 i 1   i  
2
k1 

k2 h
2

1
k 1  f t i ,  i 
k 2  f t i  h ,  i  k 1 h 
Dla i = 0, t0 = 0, Θ0 = Θ(0) = 1200:
k 1  f t 0 ,  o   f 0 ,1200    2 . 2067  10
 12
1200
4
 81  10
8
   4 .5579
k 2  f t 0  h ,  0  k 1 h   f 0  240 ,1200    4 . 5579 240   f  240 ,106 . 09 
  2 . 2067  10
1
1   0  
2
k1 
12
106 .09
4
 81  10
8
  0 .017595
1

1
k 2  h  1200     4 . 5579   0 . 017595
2
2

2
1

  240

 1200    2 . 2702 240  655 . 16 K
45
Metoda Rungego - Kutty – Przykład c.d.
Dla i = 1, t1 = t0 + h= 240, Θ1 = 655,16:
k1  f t1 ,  1   f 240 , 655 . 16    2 . 2067  10
 12
655 .16
4
 81  10
8
   0 .38869
k 2  f t1  h ,  1  k 1 h   f  240  240 , 655 . 16    0 . 38869 240   f  480 ,561 . 87 
  2 . 2067  10
1
 2  1  
2
k1 
12
561 .87
4
 81  10
8
   0 .20206

1
k 2  h  655 . 16     0 . 38869
2

2
1
1

2

    0 . 20206   240
 655 . 16    0 . 29538 240  584 . 27
Po wykonaniu dwóch iteracji metody Heun’a
otrzymujemy że temperatura kuli po 480 s
wyniesie 110.32K
46
Metoda Rungego - Kutty – Przykład c.d.
Dobór odpowiedniego kroku h jest kluczowy dla
odpowiedniej dokładności rozwiązania stosując
metodę Heun’a
temperatura Θ(K)
dokładne rozwiązanie
czas t(s)
rozmiar
kroku h
  480
480
240
120
60
30
-393.87
584.27
651.35
649.91
648.21

| t | %
160.82
9.78
0.58
0.36
0.10
47
Metoda Rungego - Kutty – Przykład c.d.
Porównanie dotychczas przedstawionych metod:
temperatura Θ(K)
Dokładna wartość rozwiązania
obliczona analitycznie wynosi:
 ( 480 )  647 . 57 K
dokładne rozwiązanie
czas t(s)
Rozmiar
kroku h
480
240
120
60
30
Θ(480)
Euler
Heun
-987.84 -393.87
110.32 584.27
546.77 651.35
614.97 649.91
632.77 648.21
Midpoint Ralston
1208.4
976.87
690.20
654.85
649.02
449.78
690.01
667.71
652.25
648.61
48
Rozwiązywanie równań różniczkowych
zwyczajnych – metoda Rungego - Kutty
Dla metody Rungego-Kutty czwartego rzędu wzór:
y i  1  y i  f  x i , y i h 
1
f
'
2!
1
 x i , y i h 2

k 1  2 k 2
 2 k 3  k 4 h
3!
f
''
 x i , y i h 3

1
4!
f
'''
 x i , y i h 4
można zapisać jako:
y i 1  y i 
1
6
gdzie najczęściej przyjmuje się że:
k1  f  x i , y i 
1
1


k 2  f  xi  h , y i  k1h 
2
2


1
1


k 3  f  xi  h , yi  k 2 h 
2
2


k 4  f  xi  h , yi  k 3h 
49
Metoda Rungego - Kutty - Przykład
Kula nagrzana do temperatury 1200 K i schładza się do
temperatury 300K. Zakładając że w procesie
schładzania ciepło jest oddawane do otoczenia
jedynie przez radiację to temperatura kuli opisana
jest równaniem:
d
dt
  2 . 2067  10
 12

4

 81  10 ,
8
 0   1200 K
Jaka jest temperatura kuli po 480 sekundach jej
schładzania?
f t ,     2 . 2067  10
 12

4
 81  10
8

50
Metoda Rungego - Kutty – Przykład c.d.
Dla metody Rungego - Kutty czwartego rzędu:
 i 1   i 
1
6
k 1
 2 k 2  2 k 3  k 4 h
Dla i = 0, t0 = 0, Θ0 = 1200:
k 1  f t 0 ,  0   f 0 ,1200    2 . 2067  10
 12
1200
4
 81  10
8
   4 .5579
1
1
1
1




k 2  f  t 0  h ,  0  k 1 h   f  0   240 ,1200    4 . 5579   240 
2
2
2
2




 f 120 , 653 . 05    2 . 2067  10
 12
653 .05
4
 81  10
8
   0 .38347
51
Metoda Rungego - Kutty – Przykład c.d.
1
1
1
1




k 3  f  t 0  h ,  0  k 2 h   f  0   240 ,1200    0 . 38347   240 
2
2
2
2




 f 120 ,1154 . 0    2 . 2067  10
 12
1154 . 0
4
 81  10
8
   3 . 8954
k 4  f t 0  h ,  0  k 3 h   f 0  240 ,1200    3 . 894   240
 f  240 , 265 . 10    2 . 2067  10
1   0 
1
6
 12
265 . 10
4
 81  10
8

  0 . 0069750
( k1  2 k 2  2 k 3  k 4 ) h
 1200 
1
6
  4 . 5579
 2   0 . 38347
  2   3 . 8954   0 . 069750 240
 1200    2 . 1848 240  675 . 65
52
Metoda Rungego - Kutty – Przykład c.d.
Dla i = 1, t1 = 240, Θ1 = 675,65:
k 1  f t1 ,  1   f  240 , 675 . 65    2 . 2067  10
 12
675 . 65
4
 81  10
8
   0 . 44199
1
1
1
1




k 2  f  t1  h ,  1  k 1 h   f  240   240 , 675 . 65    0 . 44199 240 
2
2
2
2




 f 360 , 622 . 61    2 . 2067  10
 12
622 . 61
4
 81  10
8
   0 . 31372
1
1
1
1




k 3  f  t1  h ,  1  k 2 h   f  240   240 , 675 . 65    0 . 31372   240 
2
2
2
2




 f 360 , 638    2 . 2067  10
 12
638 . 00
4
 81  10
8
   0 . 34775
k 4  f t1  h ,  1  k 3 h   f  240  240 , 675 . 65    0 . 34775   240
 f  480 , 592 . 19   2 . 2067  10
 2  1 
1
6
 12
592 . 19
4
 81  10
8

   0 . 25351
( k1  2 k 2  2 k 3  k 4 ) h
 675 . 65 
1
6
 675 . 65 
1
6
  0 . 44199
 2   0 . 31372   2   0 . 34775     0 . 25351
  2 . 0184   240
  240
 594 . 91 K
53
Metoda Rungego - Kutty – Przykład c.d.
temperatura Θ(K)
Dobór odpowiedniego kroku h jest kluczowy dla
odpowiedniej dokładności rozwiązania stosując
metodę Rungego - Kutty czwartego rzędu.
dokładne rozwiązanie
czas t(s)
rozmiar
kroku h
480
240
120
60
30
  480

| t | %
-90.278
113.94
594.91
8.1319
646.16
0.21807
647.54 0.0051926
647.57 0.00013419
54
Metoda Rungego – Kutty dla równań
różniczkowych wyższych rzędów
Dla równania różniczkowego zwyczajnego wyższego
rzędu albo dla równania cząstkowego
n
d y
an
dx
n
 a n 1
d
n 1
dx
y
n 1
   a1
dy
dx
 a o y  f x 
można dokonać podstawienia:
y  z1
dy
dz 1

dx
dx
2
d y
dx
2

 z2
dz 2
dx
 z3

d
n 1
dx
n 1
n
d y
dx
y
n

dz n 1
dx
 zn
n 1

1 
d y
dy
1
  a n 1







a

a
y

f
x

1
0
n 1
 a   a n  1 z n   a 1 z 2  a 0 z 1  f  x 
dx
a n 
dx
dx

n
dz n
55
Metoda Rungego – Kutty dla równań
różniczkowych wyższych rzędów
Otrzymane równania tworzą w efekcie układ n równań:
dz 1
dx
dz 2
dx
 z 2  f 1  z1 , z 2 ,  , x 
 z 3  f 2  z1 , z 2 ,  , x 

dz n
dx

1
an
  a n 1 z n   a 1 z 2  a 0 z 1  f  x 
Każde z n równań może być rozwiązane z użyciem
opisanych wcześniej metod rozwiązywania układów
równań różniczkowych zwyczajnych stopnia
pierwszego
56
Przykład
2
d y
Rozwiąż równanie:
dt
2
2
dy
dt
 y  e , y 0   1,
t
dy
dt
0   2
oraz oblicz y(0.75)
Podstawiając:
dy
 z
dt
Po podstawieniu równanie przybiera postać:
dz
 2z  y  e
t
dt
dz
e
t
 2z  y
dt
57
Przykład c.d.
W efekcie należy rozwiązać następujący układ równań:
dy
dt
dz
dt
 z  f 1 t,y,z , y (0)  1
e
t
 2 z  y  f 2 t , y , z , z (0)  2
Stosując metodę Eulera:
y i 1  y i  f 1 t i , y i , z i h
z i 1  z i  f 2 t i , y i , z i h
58
Przykład c.d.
Krok 1: i  0 , t 0  0 , y 0  1, z 0  2
y i 1  y i  f 1 t i , y i , z i h
y 1  y 0  f 1  t 0 , y 0 , z 0 h
 1  f 1 0 ,1, 2 0 . 25 
 1  2  0 . 25 
 1 .5
z i 1  z i  f 2 t i , y i , z i h
z 1  z 0  f 2 t 0 , y 0 , z 0 h
 2  f 2 0 ,1, 2 0 . 25 

 2 e
0

 2  2   1 0 . 25 
1
y 1  y  0 . 25   1 . 5
z 1  z 0 . 25  
dy
dt
0 . 25   1
t  t 1  t 0  h  0  0 . 25  0 . 25
59
Przykład c.d.
Krok 2:
i  1, t1  0 . 25 , y 1  1 . 5 , z1  1
y i 1  y i  f 1 t i , y i , z i h
y 2  y 1  f 1  t 1 , y 1 , z 1 h
 1 . 5  f 1 0 . 25 ,1 . 5 ,1 0 . 25 
 1 . 5  1 0 . 25 
 1 . 75
z i 1  z i  f 2 t i , y i , z i h
z 2  z 1  f 2  t 1 , y 1 , z 1 h
 1  f 2 0 . 25 ,1 . 5 ,1 0 . 25 
 1  e
 0 . 25
 2 1   1 . 5 0 . 25 
 1    2 . 7211 0 . 25 
 0.31970
y 2  y 0 . 5   1 . 75
z 2  z 0 . 5   0 . 3197
t  t 2  t 1  h  0 . 25  0 . 25  0 . 50
60
Przykład c.d.
Krok 3: i  2 , t 2  0 . 5 , y 2  1 . 75 , z 2  0 . 31970
y i 1  y i  f 1 t i , y i , z i h
z 3  z 2  f 2  t 2 , y 2 , z 2 h
y 3  y 2  f 1  t 2 , y 2 , z 2 h
 1 . 75  f 1 0 . 50 ,1 . 75 , 0 . 31970
 1 . 75  0 . 31970
0 . 25 
 1 . 8299
z i 1  z i  f 2 t i , y i , z i h
0 . 25 
 0 . 31972  f 2 0 . 50 ,1 . 75 , 0 . 31970
 0 . 31972  e
 0 . 50
 2 0 . 31970
 0 . 31972    1 . 7829
0 . 25 
  1 . 75 0 . 25 
0 . 25 
  0 . 1260
y 3  y 0 . 75   1 . 8299
z 3  z 0 . 75    0 . 12601
t  t 3  t 2  h  0 . 5  0 . 25  0 . 75
61
Przykład c.d.
Otrzymane rozwiązanie to:
y 0 . 75   y 3  1 . 8299
Wartość dokładna to:
y 0 . 75  exact  1 . 668
Błąd względny otrzymanego rozwiązania wynosi:
t 
1 . 668  1 . 8299
 100
1 . 668
62
Warunki początkowe i brzegowe
Zależność przemieszczenia v belki od długości x
oraz obciążenia q:
υ
υ
q
q
x
x
L
L
d 
2
d 
2
dx
2

q
2 EI
L  x 
2
Aby rozwiązać równanie potrzebne są dwa
warunki początkowe w punkcie x = 0
 0   0 ;
d
dx
0   0
dx
2

qx
2 EI
x  L 
Aby rozwiązać równanie potrzebne
są dwa warunki brzegowe punkach
x = 0 oraz x = L
 0   0 ;   L   0
63
Metoda różnicowa dla zagadnień brzegowych
równań różniczkowych zwyczajnych.
Zastosowanie tej metody zostanie omówione na
przykładzie równań różniczkowych drugiego rzędu
które mają narzucone warunki brzegowe.
2
d y
dx
2
 f ( x , y , y ' ), a  x  b
Warunki brzegowe:
y (a )  ya
y (b )  y b
64
Metoda różnicowa dla zagadnień brzegowych
równań różniczkowych zwyczajnych.
y
q
T
T
x
L
Odkształcenie belki wzdłuż osi Y opisuje równanie:
2
d y
dx
2

Ty
EI

qx ( L  x )
2 EI
Oblicz wartość odkształcenia belki w punkcie x = 50:
T  7200 lbs ; q  5400 lbs / in ; L  75 in ; E  30 Msi ; I  120 in
4
65
Metoda różnicowa dla zagadnień brzegowych
równań różniczkowych zwyczajnych.
2
d y
dx

2
7200 y
( 30  10 )( 120 )
6

( 5400 ) x ( 75  x )
2 ( 30  10 )( 120 )
6
2
d y
dx
Aproksymujemy
 2  10
2
2
d y
dx
2
i 1
2
d y
dx
2

6
y  7 . 5  10
7
x ( 75  x )
w punkcie i:
i
i1
y i  1  2 y i  y i 1
(x)
2
66
Metoda różnicowa dla zagadnień brzegowych
równań różniczkowych zwyczajnych.
Po podstawieniu wartości:
y i 1  2 y i  y i 1
(x)
2
 2  10
6
y i  7 . 5  10
7
x i ( 75  x i )
Ponieważ Δx = 25 będą rozpatrywane 4 węzły:
i 1
i2
i3
i4
x0
x  25
x  50
x  75
x0  0
x1  x 0   x  0  25  25
x 2  x1   x  25  25  50
x 3  x 2   x  50  25  75
67
Metoda różnicowa dla zagadnień brzegowych
równań różniczkowych zwyczajnych.
Węzeł pierwszy (i = 1):
x0
y1  0
Węzeł drugi (i = 2):
y 3  2 y 2  y1
( 25 )
2
 2  10
6
y 2  7 . 5  10
7
x 2 ( 75  x 2 )
0 . 0016 y1  0 . 003202 y 2  0 . 0016 y 3  7 . 5  10
7
( 25 )( 75  25 )
0 . 0016 y1  0 . 003202 y 2  0 . 0016 y 3  9 . 375  10
4
68
Metoda różnicowa dla zagadnień brzegowych
równań różniczkowych zwyczajnych.
Węzeł trzeci (i = 3):
y4  2 y3  y2
( 25 )
2
 2  10
6
y 3  7 . 5  10
7
x 3 ( 75  x 3 )
0 . 0016 y 2  0 . 003202 y 3  0 . 0016 y 3  7 . 5  10
7
( 50 )( 75  50 )
0 . 0016 y 2  0 . 003202 y 3  0 . 0016 y 3  9 . 375  10
4
Węzeł czwarty (i = 4):
x  75
y4  0
69
Metoda różnicowa dla zagadnień brzegowych
równań różniczkowych zwyczajnych.
Otrzymane równania dla wszystkich 4 węzłów można
zapisać jako:
 1

0 . 0016

 0

 0
0
0
 0 . 003202
0 . 0016
0 . 0016
 0 . 003202
0
0

  y1  0
  
4 
y2
0
9
.
375

10

   
0 . 0016   y 3   9 . 375  10  4 

  
1   y 4   0

0
Rozwiązanie powyższego kładu równań daje:
 y1  0
  
y
 0 . 5852
 2  
 y 3    0 . 5852
  
 y 4  0






70
Metoda różnicowa dla zagadnień brzegowych
równań różniczkowych zwyczajnych.
Ostatecznie wartość przemieszczenia y w punkcie
x = 50 wynosi:
y ( 50 )  y ( x 2 )  y 2   0 . 5852 "
Rozwiązanie analityczne daje:
y  0 . 375 x  28 . 125 x  3 . 75  10  1 . 775656266  10 e
2
5
5
0 . 0014142 x
 1 . 974343774  10 e
5
 0 . 0014142 x
y ( 50 )   0 . 5320
Błąd rozwiązania wyznaczonego metodą różnicową:
t 
 0 . 5320  (  0 . 5852 )
 0 . 5320
 100 %
 t  10 %
71
Mikro- i nanobelki w sensorach
72
Zagadnienie problemowe – równanie
przewodnictwa cieplnego
Równanie przewodnictwa cieplnego zwane jest także jako
równanie dyfuzji.
W celu wyprowadzenia tego równania rozważamy podobszar
V obszaru Ω, o gładkim brzegu ∂V.
Niech F oznacza gęstość strumienia przepływu w obszarze
Ω, wówczas tempo w jakim zmienia się ilość substancji
wypełniającej V jest równe całkowitemu przepływowi
substancji przez brzeg ∂V:
d
dt
u
(
t
,
x
)
dx



V

V
F  dS    divFdx
V
η – wektor normalny do ∂V, dS - miara na powierzchni ∂V
73
Równanie przewodnictwa cieplnego –
rozwiązanie podstawowe
Równanie ciepła ma bardziej skomplikowaną strukturę do
równania Laplace’a – poszukiwanie rozwiązań w postaci
tzw. funkcji samopodobnej tzn.:
 x 
n
u (t , x )   v     ,   R , R  R
t
t 
1
Podstawiając do równania cieplnego otrzymujemy:
t
  1
yt
v( y)  t

  1
y  v( y)  t
  2 
v( y)  0
x
74
Równanie przewodnictwa cieplnego –
rozwiązanie podstawowe
Jeżeli β=1/2 nasze równanie ulega uproszczeniu:
v( y) 
1
y  v( y)  v( y)  0
2
Stosując kolejne uproszczenie i operator Laplace’a:
v ( y )  w ( r ) gdzie r  y oraz w : R  R
 w(r ) 
1
2
r  w (r )  w (r ) 
'
''
n 1

?
w (r )  0
'
r
75
Równanie przewodnictwa cieplnego –
rozwiązanie podstawowe
Jeżeli przyjmiemy, że: α=n/2 to otrzymamy równanie:
n
2
w(r ) 
1
r  w (r )  w (r ) 
'
''
n 1
2
w (r )  0
'
r
Stosując pewne techniki mnożenia i dzielenia rn-1
otrzymujemy:
 r 

w ( r )  C 2 exp  

4


2
76
Równanie przewodnictwa cieplnego –
rozwiązanie podstawowe
I ostatecznie uwzględniając wszystkie założenia:
v ( y )  w ( y ), y  t

x,   n 2,   1 2
Otrzymujemy funkcję, która jest rozwiązaniem równania
ciepła:
2

x 
C
n


exp 
, t  0, x  R
n/2
 4t 
t


77
Równanie przewodnictwa cieplnego –
rozwiązanie podstawowe
Rozwiązaniem podstawowym operatora
przewodnictwa cieplnego nazywamy funkcję:
2

 x 
1
n



exp 
dla t  0 , x  R ,
n 2
 4t 
E ( x , t )   ( 4 t )



n
dla t  0 , x  R
 0
78
Rozkład temperatury w czujnikach
79
Zagadnienie początkowe dla równania struny.
Wzór d’Alamberta
Badanie równania falowego zaczniemy od
przypadku jednowymiarowego czyli od tzw.
Równania struny.
Skoncentrujemy się na równaniu struny
nieograniczonej.
Naszym celem jest rozwiązanie zagadnienia:
 u tt  c 2 u xx  0 dla t  0 , x  R ,

dla x  R ,
 u (0, x )  g ( x )
 u (0, x )  h ( x )
dla
x

R
,
t

80
Zagadnienie początkowe dla równania struny.
Wzór d’Alamberta
Rozwiązanie ogólne równania wyraża się
wzorem:
u ( t , x )  F ( x  ct )  G ( x  ct )
Gdzie F i G są funkcjami klasy C2(R). Korzystając
z warunków początkowym dostajemy:
 u (0, x )  F ( x )  G ( x )  g ( x )
 '
'
'
 u t ( 0 , t )  cF ( x )  cG ( x )  h ( x )
81
Zagadnienie początkowe dla równania struny.
Wzór d’Alamberta
Całkując drugie równanie mamy:
F ( x)  G ( x) 
1

c
x
x0
h ( r ) dr  F ( x 0 )  G ( x 0 )
Zatem rozwiązaniem powyższego układu równań
jest para funkcji:
F ( x) 
1

2
G ( x) 
1
2
1
2c

1
2c
x
 h ( r ) dr

2
x0
x
 h ( r ) dr
x0
F ( x0 )  G ( x0 )

F ( x0 )  G ( x0 )
2
82
Zagadnienie początkowe dla równania struny.
Wzór d’Alamberta
Stąd rozwiązanie problemu struny wyraża się
wzorem d’Alamberta:
u (t , x ) 
1
2
( g ( x  ct )  g ( x  ct )) 
1
2c
x  ct
h
(
r
)
dr

x  ct
Jeśli gєC2(R), hєC1(R), to rozwiązanie zagadnienia
struny wyraża się wzorem d’Alamberta
Zadanie domowe: wymuszone drgania struny
83
Rezonans struny
84
Określenie stabilności wg Łapunowa
Rozważmy układ macierzowy równań
różniczkowych w postaci:

x  Ax
Pytanie: Jaki punkt jest stabilny dla układu
liniowego?
85
Określenie stabilności wg Łapunowa
Rozważmy:
x '  Tx
A'  T
1
AT  diag ( 1 ,...,  n )
λ – wartości własne macierzy
det( A   I )  0
86
Określenie stabilności wg Łapunowa
Otrzymujemy:
xi ' ~ e
xT
 it   i  0
1
 0
x' 0
x – jest to punkt asymptotycznie stabilny
87
Określenie stabilności wg Łapunowa
Jeżeli wartości własne macierzy A są mniejsze od
zera wówczas możemy powiedzieć, że:
  , r  0  (t ,  )  r  e
 t
t 0
x – jest to punktem stabilnym wg Łapunowa i
asymptotycznie stabilnym
88
Określenie stabilności wg Łapunowa
Punkt nazywamy stabilnym wg Łapunowa jeżeli
spełnione są warunki (3):
   0   a  p  y (t ,  )  t
   0    p   a     (t ,  )  a    t
a  p .stab .asymptotyc znej
  
  a   oraz lim  ( t ,  )  a  0
t  
89
Określenie stabilności wg Łapunowa
Sens definicji Łapunowa ilustruje rysunek:
x ( t )  f ( x ( t ), t )
90
Określenie stabilności wg Łapunowa
Lokalna asymptotyczna stabilność wg Łapanowa
x ( t )  f ( x ( t ), t )
91
Stabilność rozwiązań równań różniczkowych
Punkt x0 nazywamy punktem stacjonarnym (położeniem równowagi)
równania:
x '  f ( x )  f ( x0 )  0
Stabilność w sensie Łapunowa – jeśli startując z warunku
początkowego x0 blisko rozwiązania stacjonarnego, pozostajemy w
pobliżu tego rozwiązania wraz z upływem czasu.
Asymptotycznie stabilne – jeśli jest stabilne i dla warunku
początkowego x0 dostatecznie blisko x=const., rozwiązanie x(t) z tym
warunkiem początkowym zbiega do x przy t->∞.
Niestabilne r.r. – jeśli znajdzie się taki punkt początkowy dowolnie
blisko x=const., dla którego rozwiązanie ucieka wraz z upływem
czasu.
92
Stabilność rozwiązań równań różniczkowych
Twierdzenie Hartmana-Grobmana
x '  Ax ,  ( A )  spectrum
Jeżeli f(p)=0 i wśród wartości własnych Df(p) nie ma wartości własnych
czysto urojonych to równanie nieliniowe x’=f(x) i liniowe x’=df(p) są
topologicznie sprzężone w otoczeniu p.
Tw. Łapunowa
Re(  ( A ))  0  x  0
Asymptotycznie stabilny pkt. równowagi
Re(  ( A ))  0  x  0
Stabilny punkt równowagi
    ( A )  Re(  )  0  x  0
Niestabilny punkt
93
Stabilność rozwiązań równań różniczkowych przykład
Rozważmy układ równań:
 y'  y

2
z' z  y
Szukamy punktów równowagi:
 y  0
y  0




0
,
0
 2

2
z  0
z  y  0
94
Stabilność rozwiązań równań różniczkowych przykład
Macierz linearyzacji:
  f1

 dy
 f 2

 dy
 f1 

 1
dz 
 
f 2 
2y

dz 
0

1
Obliczamy macierz linearyzacji w podanym punkcie (0,0):
1

 0
0

1
Wartości własne to: -1;1. Nie ma wartości czysto urojonych więc
punkt (0,0) jest niestabilny.
95
Metoda różnic skończonych
Metoda różnic skończonych opiera się na zastąpieniu
pochodnych cząstkowych w punktach (xi,yi) ich przybliżeniami
numerycznymi. Otrzymujemy odpowiednio dla zmiennej x i y:
 u
2
x
2
 u
( xi , y i ) 
2
y
2
( xi , y i ) 
u ( x i  1 , y j )  2 u ( x i , y j )  u ( x i 1 , y j )
h
2
u ( x i , y j  1 )  2 u ( x i , y j )  u ( x i , y j 1 )
k
2
h  u
2

12  x
4
h  u
2

4
( i , y j )
4
12  y
4
( x i , j )
 i  ( x i 1 , x i  1 ) oraz  i  ( y j 1 , y j  1 )
96
Metoda różnic skończonych dla równania
Poissona
Podstawienie FEM do równania Poissona otrzymujemy:
u ( x i  1 , y j )  2 u ( x i , y j )  u ( x i 1 , y j )
h
h  u
2
 f ( x, y ) 
2
k
4
12  x
4

u ( x i , y j  1 )  2 u ( x i , y j )  u ( x i , y j 1 )
( x i , j ) 
dla i  1, 2 ,..., n  1 oraz
k
2
 u
2
4
12  y
4
( x i , j )
j  1, 2 ,..., m  1
Warunki brzegowe:
u ( x 0 , y j )  g ( x 0 , y j ) oraz u ( x n , y j )  g ( x i , y j ) dla j  0 ,1,..., m
u ( x i , y m )  g ( x i , y m ) oraz u ( x i , y 0 )  g ( x i , y 0 ) dla i  0 ,1,..., n
97

Metoda różnic skończonych dla równania
Poissona
Pomijając reszty otrzymujemy układ (n-1) x (m-1) równań liniowych
z niewiadomymi, które są przybliżeniami u(xi,yj).
Układ równań możemy rozwiązać metodami bezpośrednimi bądź
iteracyjnymi.
W celu wyznaczenia
przybliżenia w punkcie
(xi,yj), potrzebne są
wartości przybliżenia
rozwiązania w czterech
sąsiednich punktach:
98
Metoda różnic skończonych przykład
Wyznaczyć rozkład temperatury w stanie ustalonym dla cienkiej
kwadratowej metalowej płytki o wymiarach 0,5m na 0,5m.
Na brzegu płytki znajdują się źródła ciepła utrzymujące temperaturę
na poziomie 0oC dla boku dolnego i prawego, natomiast temperatura
boku górnego i lewego zmienia się liniowo od 0oC do 100oC.
Problem rozwiązać układając układ równań liniowych (postać
macierzowa) dla wewnętrznych węzłów siatki 5 x 5 – układ równań
rozwiązać metodą Gaussa-Siedla.
99
Metoda różnic skończonych przykład
100
Metoda różnic skończonych przykład
Problem ten opisuje równanie Laplace’a
2
 T ( x, y ) 
2
d T
dx
2
2
( x, y ) 
d T
dy
2
( x, y )  0
Z warunkami brzegowymi:
T ( 0 , y )  0 [ C ] T ( x ,0 )  0 [ C ]


T ( 0 , 5 , y )  200 y [ C ] T ( x , 0 ,5 )  200 x [ C ]


101
Metoda różnic skończonych przykład
Postać macierzowa układu: [zadanie domowe – obliczyć temperatury]
102
Metoda różnic skończonych – równania
paraboliczne
Przypomnijmy równanie paraboliczne:
du
dt
2
( x, t )  
2
d u
dx
2
( x, t )  0
dla 0  x  l oraz t  0
Z warunkami brzegowymi i początkowymi:
u ( 0 , t )  u (l , t )  0
dla t  0
u ( x ,0 )  f ( x )
0 xl
103
Metoda różnic skończonych – równania
paraboliczne
Dla danego m definiujemy krok h=(b-a)/m. Ustalamy wartość kroku
czasowego k. Stąd węzły siatki (xi,tj):
x i  ih
 u
2
x
2
 u
( xi , y i ) 
2
t
2
( x i ,t i ) 
dla i  0 ,1,..., m oraz
t j  jk dla
u ( x i  1 , y j )  2 u ( x i , y j )  u ( x i 1 , y j )
h
2
u ( x i , t j  1 )  2 u ( x i , t j )  u ( x i , t j 1 )
k
2
j  0 ,1,... m
2

4
12  x
h  u
2

h  u
4
( i , y j )
4
12  t
4
( x i , j )
 i  ( x i 1 , x i  1 ) oraz  i  ( t j 1 , t j  1 )
104
Metoda różnic skończonych – równania
paraboliczne
Otrzymujemy:
u ( x i , t j 1 )  u ( x i , t j )

2
u ( x i  1 , t j )  2 u ( x i , t j )  u ( x i 1 , t j )
k
h
2
Po uwzględnieniu warunku brzegowego: u ( 0 , t j )  u ( l , t j )  0
w i , j 1  w i , j
k

2
w i  1, j  2 w i , j  w i 1, j
h
2
0
dla j  1,.., m
0
w i , j  1  (1  2  ) w i , j   ( w i  1 , j  w i 1 , j )
2
 
 , i  1, 2 ,..., m  l ,
 h 
  k
j  1, 2 ,... m
105
Metoda różnic skończonych – równania
paraboliczne
Schemat jawny
Warunek
zbieżności
schematu
jawnego

2
k
h
2

1
2
106
Metoda różnic skończonych – równania
paraboliczne
w i , j  w i , j 1
Schemat
niejawny:
k

2
w i  1, j  2 w i , j  w i 1, j
h
2
0
w i , j 1  (1  2  ) w i , j   ( w i  1 , j  w i 1 , j )
2
 
 , i  1, 2 ,..., m  l ,
 h 
  k
j  1, 2 ,... m
Schemat niejawny
jest zawsze zbieżny,
niezależnie od
wielkości kroku
całkowania
107