Transcript Document

Metody matematyczne w
Inżynierii Chemicznej
Rozwiązywanie równań
różniczkowych
Równanie różniczkowe rzędu n
Wzór ogólny

f x, y, y, y,..., y
n 
 0
Cel rozwiązania równania
różniczkowego
Matematyk: rozwiązanie ogólne w
postaci funkcji
Inżynier: rozwiązanie szczególne,
określone przez warunki początkowe.
Interesujące są wartości funkcji dla
kolejnych zmiennych niezależnych, czyli
zbiór par (x1, y1), (x2,y1),...,(xn, yn) kiedy
dana jest funkcja f(x, y, y', y",..,y(n) )=0
Warunki początkowe
Zagadnienie początkowe
Zagadnienie brzegowe
Warunki początkowe
Zagadnienie początkowe – wszystkie
równania warunków początkowych
podane są dla tej samej zmiennej
niezależnej
Np. dla równania:
d2y
dy
 a  by  0
2
dt
dt
Warunki początkowe:
t0
 dy 
  0
 dt t 0
 yt 0  0
Warunki początkowe
Zagadnienie brzegowe –równania
warunków początkowych podane są dla co
najmniej dwóch wartości zmiennej
niezależnej
Np. dla równania:
Warunki brzegowe:
d2y
dy
 a  by  0
2
dx
dx
x0
 dy 
  0
 dx  x 0
x0
 y x0  0
lub
i
i
x  10
 y x10  1
x  10
 y x10  1
Równania wyższych rzędów
Przekształca się do układów równań
rzędu pierwszego
Np. w równaniu:
podstawmy:
stąd:
d2y
dy
 a  by  0
2
dt
dt
dy
z
dt
d 2 y dz

2
dt
dt
Równania wyższych rzędów
Otrzymujemy układ równań pierwszego rzędu:
 dz
 dt  az  by  0

 dy  z
 dt
Równania wyższych rzędów
Dla równanie trzeciego rzędu
d3y
d2y
dy
 a 2  b  cy  0
3
dt
dt
dt
dy
z
dt
d 2z
dz
 a  bz  cy  0
2
dt
dt
Równania wyższych rzędów
Równanie trzeciego rzędu przechodzi w układ 3 równań
pierwszego rzędu
dy
z
dt
dz
x
dt
dx
 ax  bz  cy  0
dt
Metody rozwiązywania r.r.
1.
2.
Metody wielokrokowe: yi+1 oblicza się na
podstawie znanych yi, yi-1, yi-2,.., yi-p. Do
wyliczenia punktu (xi+1, yi+1) wymagana jest
znajomość p+1 punktów obliczonych
wcześniej
Metody klasy Rungego-Kutty: yi+1 oblicza się
na podstawie yi i pewnych wartości
pośrednich F(xi+a, yi+b), gdzie
a należy do przedziału <0, h>
b oblicza się wg algorytmu danej metody
Metoda Eulera
x  x0
dy
 f  x, y  z warunkami początkowymi
dx
y y

 f  x, y 
x h
 yxx
0
 yi 1  yi   f x , y   Oh 
i
i
h
yi1  yi  hf xi , yi   Oh
x0 , y0
y1 , x1  x0  h
2
 y0
II rząd

y2 , x2  x1  h
Metoda Eulera dla równań 2-go rzędu
d2y
dy 

 f  x, y, 
2
dx
dx 

 dy
 dx  z

 dz  f  x, y , z 
 dx
z warunkami
początkowymi
x  x0
 yxx0  y0
 dy 
 y0
 
 dx  x x0
x  x0
z warunkami
początkowymi
 yxx
0
 y0
z xx
 y0
0
Metoda Eulera dla równań 2-go rzędu
 yi 1  yi  hz i

 zi 1  zi  hf  xi , yi , z1 
x0
x1  x0  h x2  x1h
y0
y0  z0
y1
y2
z1
z2
Metody wielokrokowe
Typy
1.
2.
Wykorzystujące wzory na wartość
pochodnej w punkcie
Wykorzystujące metody całkowania
numerycznego
Zasada metod wielokrokowych
xi-p
xi-3
xi-2
xi-1 xi
xi+1
Metody wielokrokowe typ 1.
Pochodną w równaniu:
y  F x, y 
Podstawia się odpowiednim wyrażeniem.
Jeżeli zastosować najprostszy wzór na pochodną:
yi 1  yi
yi 
 Oh 
h
Po podstawieniu:
yi 1  yi
F  xi , yi  
 Oh 
h
i przekształceniu:
2



yi 1  yi  hF xi , yi  O h 
m. Eulera
Metody wielokrokowe typ 1.
Jeżeli zastosować dokładniejszy wzór na pochodną:
 
1
 yi 1  yi 1   O h 2
yi 
2h
Po podstawieniu:
 
1
2
 yi 1  yi 1   O h
F xi , yi  
2h
i przekształceniu:
 
yi 1  yi 1  2hFxi , yi   O h3
Metody wielokrokowe typ 1.
Jeszcze większą dokładność otrzyma się stosując wzór:
 
1
 yi 2  6 yi 1  3 yi  2 yi 1   O h3
yi 
6h
Po podstawieniu:
 
1
 yi 2  6 yi 1  3 yi  2 yi 1   O h3  F xi , yi 
6h
i przekształceniu:
 
3
1
4
yi 1   yi  3 yi 1  yi  2  3hF xi , yi   O h
2
2
Metody wielokrokowe typ 1.
Podsumowanie
 
 2hFx , y   Oh 
yi 1  yi  0 yi 1  hFxi , yi   O h
2
yi 1  0 yi  yi1
i
i
3
1
yi 1   yi  3 yi 1  yi  2  3hF xi , yi   Oh 4 
2
2
3
yi 1  a0 yi  a1 yi 1  a2 yi 2  ... a p yi  p  b1F xi , yi 
Ogólny wzór metod wielokrokowych jawnych
Metody wielokrokowe typ 2.
Opierając się na operacji całkowania równania:
y  F x, y 
W granicach przedziału <i-p, i+1>
i 1
i 1
i p
i p
 ydx   F x, y dx
Lewa strona jest dokładnie równa różnicy wartości funkcji
między punktami i-p, i+1:
i 1
yi 1  yi  p 
 F x, y dx
i p
Metody wielokrokowe typ 2.
Prawą stronę oblicza się całkując numerycznie jedną
z metod.
Jeżeli zastosować metodę prostokątów to p = 0
i 1
yi 1  yi 
yi 1  yi 
xi  h
 F x, y dx
i
 
2




F
x
,
y
dx

hF
x
,
y

O
h
i
i

xi
 
yi 1  yi  hFxi , yi   O h2
m. Eulera!!
Metody wielokrokowe typ 2.
Jeżeli zastosować metodę trapezów to p = 0
xi  h
 
h
3
yi 1  yi   F  x, y dx  F xi , yi   F xi 1 , yi 1   O h
2
xi
h
3
yi 1  yi  F xi , yi   F xi 1 , yi 1   O h
2
 
Równanie to jest
uwikłane ze
względu na yi+1.
Metody takie nazywane są niejawnymi
Metody wielokrokowe typ 2.
F xi , yi  to wartość pochodnej w punkcie i
Można ją oznaczyć yi , co upraszcza zapis
Ponieważ
Metodę bazującą na całkowaniu metodą trapezów
można ostatecznie zapisać
 
h
yi 1  yi   yi  yi1   O h 3
2
Rozwiązanie wymaga wykonania obliczeń iteracyjnych
Metody wielokrokowe typ 2.
1.
2.
3.
4.
Wstępne oszacowanie wartości yi+1.
Obliczenie pochodnej y'i+1=F(xi+1, yi+1)
Obliczenie yi+1 z wyprowadzonego wzoru
Porównanie oszacowanej i obliczonej
wartości yi+1 . Jeżeli różnią się o więcej niż
założona wartość e to powrót do punktu 2.
Metody oparte o wzory całkowe mają większą
dokładność niż bazujące na równaniach na
obliczenie pochodnych
Metody wielokrokowe typ.2
Stosując wzór całkowy Simpsona (p = 1)
i 1
i 1
 
h
5




i1y dx  i1F x, y dx  3  yi1  4 yi  yi1   O h
Otrzymuje się wzór niejawny
 
h
yi 1  yi 1   yi1  4 yi  yi1   O h 5
3
Metody wielokrokowe typ.2
podsumowanie
yi 1  yi 0  hyi0
p=0
h
h
yi 1  yi 0  yi0  yi1
2
2
p=0
h
4h
h
yi 1  yi 1  yi1 
yi  yi1
3
3
3
yi1  yi  p  bp1 yi p  bp yi p1  ... b0 yi1
p=1
Metody wielokrokowe wzór
ogólny
yi 1  a0 yi  a1 yi 1  a2 yi2  ... a p yi  p  b1 yi
yi1  yi  p  bp1 yi p  bp yi p1  ... b0 yi1
yi 1  a0 yi  a1 yi 1  ... a p yi  p  b0 yi1  b1 yi  .... bp yi p
Jest to ogólny wzór na metody wielokrokowe
b0 = 0 to wzór jest jawny,
b0  0 wzór jest niejawny
Metody wielokrokowe dwuetapowe
Pierwszy krok można wykonać stosując metodę
jawną opartą o wzory na pochodną. Np.:
1.
 
y *i 1  yi 1  2hFxi , yi   O h3
Pierwsze przybliżenie y*i+1 jest nazywane PROGNOZĄ
(PREDICTOR)
2. Obliczenie pochodnej y'i+1=F(xi+1 , y*i+1 )
3. Lepsze przybliżenie yi+1
 
h
yi 1  yi   yi1  yi   O h 3
2
Nazywane korektą (uściśleniem) - CORRECTOR
Metody wielokrokowe dwuetapowe
Punkty 2. i 3. powtarzane są iteracyjnie.
Warunek zakończenia obliczeń iteracyjnych można
przedstawić następująco:
yik1  yik11  h n
Ten sposób obliczeń nazywany jest
metodą predictor-corrector.
Przedstawiona metoda nosi nazwę
zmodyfikowanej metody Eulera
Metody wielokrokowe predictor-corrector
Metoda Milne'a (1)
Prognoza:
 
4h
2 yi  yi1  2 yi2   O h5
yi 1  yi 3 
3
Korekta:
 
1
yi 1  9 yi  yi  2  3h yi1  2 yi  yi1   O h 5
8
Metody wielokrokowe predictor-corrector
Metoda Milne'a (2)
Prognoza:
 
3h
7





yi 1  yi 5  11 yi  14 yi 1  26 yi  2  14 yi 3  11 yi  4   O h
10
Korekta:
 
2
yi 1  yi 3  h7 yi1  32 yi  12 yi1  32 yi 2  7 yi3   O h 7
45
Metody wielokrokowe predictor-corrector
Wzory Adamsa
Prognoza (Adamsa-Bashfortha):
 
h
yi 1  yi  55 yi  59 yi1  37 yi 2  9 yi3   O h 5
24
Korekta (Adamsa-Moultona):
 
h
yi 1  yi  9 yi1  19 yi  5 yi1  yi 2   O h 5
24
Obliczanie punktów początkowych w
metodach wielokrokowych
Stosując metody wielokrokowe mamy dany
warunek początkowy typu zagadnienie początkowe
czyli współrzędne punktu początkowego x0 i y0.
Aby wykonać obliczenia metodą o pewnej wartości p
potrzeba jeszcze p par xi, yi. Np. W pierwszej metodzie
Milne'a p = 3. Pierwszą wartość jaką możemy obliczyć jest
y4 a innym sposobem trzeba obliczyć y1, y2, y3.
Można wykorzystać rozwinięcie w szereg Taylora w
otoczeniu punktu x0.
1 2 2
1 3 3
1 4 4 4 
yi  y x  ih   y0  ihy0  i h y0  i h y0  i h y0  .......
2
6
24
Stabilność i zbieżność obliczeń
y2
y3
Y
y1
y0
h
x0
x1
x2
x3
Zbieżność oznacza, że:
lim y  x, h   Y  x 
h 0
Stabilność i zbieżność obliczeń
We wszystkich wzorach błąd jest dodatnią potęgą kroku:
Ponieważ krok jest bardzo mały h << 1
oraz
h1  h
2
h1
stąd
n
 hn
2
 
n

O h 1   O h n
 2 
Przedstawione metody spełniają warunek zbieżności.
Stabilność i zbieżność obliczeń
Definicja stabilności
Rozwiązanie numeryczne równania różniczkowego jest
stabilne, jeżeli błąd wniesiony do obliczeń przez
zaokrąglenie lub metodę zostanie w trakcie obliczeń
stłumiony lub rośnie wolniej od obliczonych wartości (błąd
względny maleje)
Stabilność i zbieżność obliczeń
Stabilność, tak jak zbieżność, zależy od kroku.
Błąd k-tego kroku obliczeń numerycznych to suma
błędów metody i zaokrąglenia:
ek  emk  ezk
cz
e 
h
k
z
emk  cmhn1
Stabilność i zbieżność obliczeń
0.25
błąd zaokrąglenia
błąd metody
błąd wypadkowy
0.2
0.15
0.1
0.05
0
0
0.002
0.004
hopt
0.006
0.008
0.01
Algorytm rozwiązywania równań różniczkowych
metodą jawną O(h3)
1. Czytaj parametry punktów startowych x0, y0, x1, y1
2. Czytaj końcową wartość xk i krok h
3. Podstaw za i wartość 2
4. Oblicz xi = x0+i*h
5. Oblicz yi = yi-2+2hF(xi-1, yi-1)
6. Zwiększ i o 1
7. Oblicz xi = x0+i*h
8. Jeżeli xi <= xk to idź do punktu 5
9. Podstaw za n wartość i-1
10. Podstaw za i wartość 0
11. Drukuj xi oraz yi
12. Zwiększ i o 1
13. Jeżeli i<=n to idź do punktu 11
14. Koniec
Algorytm rozwiązywania równań różniczkowych
metodą typu predictor-corrector
1. Czytaj parametry punktów startowych x0, y0, x1, y1
2. Czytaj końcową wartość xk oraz krok h
3. Podstaw za i wartość 1
4. Oblicz xi+1 = x0+(i+1)*h
5. Oblicz yi+1 = yi-1+2hF(xi, yi)
6. Przyjmij y*= yi+1
 

7. Oblicz yi 1  yi  h / 2 F xi 1 , y *  F
8. Jeżeli |y* – yi+1|>h3 to idź do punktu 6
9. Zwiększ i o 1
10. Oblicz xi+1 = x0+(i+1)*h
11. Jeżeli xi+1 <= xk to idź do punktu 5
xi , yi 
Algorytm rozwiązywania równań różniczkowych
metodą typu predictor-corrector
12. Podstaw za n wartość i-1
13. Podstaw za i wartość 0
14. Drukuj xi oraz yi
15. Zwiększ i o 1
16. Jeżeli i<=n to idź do punktu 14
17. Koniec