Wymiana ciepła w zbiornikach (1)
T0 = 20 C
Tpary = 250 C
W = 100 kg/min
M = 1000 kg
Cp = 2.0 kJ/kg
UA = 10 kJ/(min C)
Trzy zbiorniki połączone szeregowo są używane do wstępnego podgrzewania ropy
naftowej przed podaniem jej do kolumny destylacyjnej. Każdy zbiornik w chwili
początkowej jest napełniony 1000 kg ropy w temperaturze 20 C. Zawartość
zbiorników jest dobrze mieszana. Para nasycona o temperaturze 250 C
kondensuje w wężownicach zanurzonych w zbiornikach. Ciepło właściwe ropy
wynosi 2.0 kJ/kg. Szybkość przepływu ciepła w każdym zbiorniku dana jest
zależnością
Q =UA(Tpary -Ti )
Jak zmienia się w czasie temperatura ropy w poszczególnych zbiornikach?
Wymiana ciepła w zbiornikach (2)
Bilans energii dla pierwszego zbiornika
dT1
MC =WC (T0 -T1) +UA(Tpary -T1)
p
dt
Układ równań opisujący przepływ ciepła w zbiornikach
dT1
= [WC (T0 - T1) +UA(Tpary - T1)]/(MCp )
p
dt
dT2
= [WC (T1 - T2) +UA(Tpary -T2)]/(MCp )
p
dt
dT3
= [WC (T2 - T3) +UA(Tpary - T3)]/(MCp )
p
dt
Warunki początkowe
T1(0) = 20, T2(0) = 20, T3(0) = 20
Wymiana ciepła w zbiornikach (3)
function dTdt = zbiorniki (t, T, T0, W, M, Cp, Tpary, UA)
dTdt = [(W*Cp*(T0-T(1)) + UA*(Tpary-T(1)))/(M*Cp)
(W*Cp*(T(1)-T(2)) + UA*(Tpary-T(1)))/(M*Cp)
(W*Cp*(T(2)-T(3)) + UA*(Tpary-T(1)))/(M*Cp)];
Temperatura w zbiornikach
55
50
uchwyt funkcji
45
40
W = 100; 35
M = 1000;
30
W = 100;
25
Cp=2;
Tpary=250;
20
0 20 40 60 80 100
UA=10;
czas [min]
T0=20;
[t, T] = ode45(@(t, T) zbiorniki(t, T, T0,W, M, Cp,Tpary,UA), [0 90], ones(1,3)*T0);
o
temperatura [ C]
Manometr (1)
Jak będzie zmieniał się poziom wody w
manometrze gdy na jednym z jego końców
ciśnienie zwiększy się o 1000 Pa?
W chwili początkowej manometr znajduje
się w stanie równowagi, a ciśnienia na jego
obu końcach są takie same. Dane
gęstość wody: r = 1000 kg/m3
lepkość wody: m = 0.001239 Pas
promień rurki: R = 0.005 m
długość cieczy w rurce: L =1 m
przyspieszenie ziemskie: g = 9.81 m/s2
Warunki początkowe
2
L d h 4mL dh Dp
+ + h =
dh
2g dt 2rg
dt2 rg R2
h(0) = 0, = 0
dt
t=0
Manometr (2)
2
L d h 4mL dh Dp
+ + h =
dx1
2g dt 2rg
dt2 rg R2
= x2
dt
2
d h 8m dh 2g 1
dx2
+ + h = Dp
= cDp - ax2 - bx1
dt L rL
dt2 r R2
dt
8m 2g 1
Warunki początkowe
a = , b = , c =
L rL
r R2
x1(0) = 0, x2(0) = 0
2
d h dh
+ a + bh = cDp
dt
dt2
x1 = h
function dxdt = manometr (t, x ,a, b, c, dp)
dh
dxdt = [x(2)
x2 =
dt
c*dp - a*x(2) - b*x(1)];
Manometr (3)
Manometr
0.3
wychylenie [m]
predkosc [m/s]
0.2
ro = 1000;
0.1
mi = 0.001239;
L = 1; 0
R = 0.005;
-0.1
g = 9.81;
-0.2
0 5 10 15 20
a = 8*mi/(ro*R^2);
czas [s]
b = 2*g/L;
c = 1/(ro*L);
Dp = 1000;
[t, x]=ode45(@(t, x) manometr (t, x, a, b, c, Dp), [0 20], [0 0]);
Przekształcenie do postaci kanonicznej
n 2 n-1
ć
d y dy d y d y
= f y, , ,L, , x
dx
dxn dx2 dxn-1 ł
Ł
y = z1
dz1
dy dz1
= z2
== z2
dx
dx dx
2
dz2
d y d dy dz2
ć
= z3
= =
dx = z3 dx
dx2 dx dx
Ł ł
M
n-1 n-2
dzn-1
ć
d y d d y dzn-1
= zn
= = = zn
dx
dxn-1 dx dxn-2 dx
Łł
dzn
nn-1
ć
d y d d y dzn
= f (z1,z2,z3, K, zn,x)
==
dx
dxn dx dxn-1 dx
Łł
Zastosowanie wzoru Taylora
x(n+1) (t +qh)
óó Rn = hn+1
x (t) x(n)(t)
ó
x(t + h) = x(t) + x (t)h + h2 + ...+ hn + Rn
(n +1)!
2! n!
gdzie 0 < q <1
óó
x (t) x(n)(t)
ó
- metoda rzędu n-tego
x(t + h) x(t) + x (t)h + h2 + ...+ hn
2! n!
Błąd lokalny metody (spowodowany odrzuceniem składników powyżej n-tego):
x(n+1)(t +qh)
Rn = hn+1 gdzie 0
(n +1)!
Błędy rozwiązywania
Biorąc najprostsze przybliżenie n+1 pochodnej numerycznego r.r.:
x(n)(t + h) - x(n)(t) - błąd lokalny metody
x(n+1)(t)
- błąd lokalny zaokrąglenia
h
- błąd globalny metody
1
Rn hn[x(n)(t + h) - x(n)(t)]
- błąd globalny zaokrąglenia
(n +1)!
- błąd całkowity
błąd globalny suma błędów lokalnych
błąd całkowity suma błędów globalnych metody i zaokrąglania
Zastosowanie wzoru Taylora. Przykład.
dx
= cos t - sin x + t2, x(-1) = 3
dt
óó óóó
x (t) x (t) x(4)(t)
ó
x(t + h) x(t) + x (t)h + h2 + h3 + h4 - metoda rzędu 4-tego
2! 3! 4!
óó ó
x = -sin t - x cosx + 2t
óóó óó ó
x = -cost - x cosx + (x )2 sin x + 2
Zalety:
- prostota
óóó ó óó ó
x(4) = sin t - x cosx + 3x x sin x + (x )3 cosx
- możliwa do osiągnięcia
input M Ź 200; h Ź 0.01; t Ź -1.0; x Ź 3.0
wysoka dokładność
output 0, t, x
Wady:
for k = 1 to M do
ó - wielokrotne różniczkowanie
x Ź cost - sin x + t2
óó ó
x Ź -sin t - x cosx + 2t
funkcji f
óóó óó ó
x Ź -cost - x cosx + (x )2 sin x + 2
óóó ó óó ó
x(4) Ź sin t - x cosx + 3x x sin x + (x )3 cosx
1 1 1
- schemat Hornera
ó óó óóó
x Ź x + h(x + h(x + h(x + hx(4) )))
2 3 4
t Ź t + h
output k, t, x
- 1 < t Ł 1
end do
Metoda Eulera
Na podstawie wzoru Taylora
dx
= f (t, x)
dt
ó
x(t + h) = x(t) + x (t)h + R1
ó
x(t + h) x(t) + x (t)h
x(t + h) x(t) + f (t, x)h
wartość w następnym kroku = wartość w poprzednim kroku + nachylenie * krok
Błąd lokalny metody
óó
x (t +qh)
R1 = h2 = O(h2)
2!
Błąd globalny metody jest O(h)
Metoda Eulera. Przykład
8
dx
= -2t3 +12t2 - 20t + 8.5
analityczne
7
dt
h = 0.5
h = 0.25
6
x(0) =1
5
4
t = 0;
tf = 4;
3
h = input ('h = ');
2
x = 1;
wt = t;
1
0 1 2 3 4
wx = x;
t
while t < tf
dxdt = -2*t^3 + 12*t^2 -20*t + 8.5;
x = x + dxdt*h;
t = t + h;
wt = [wt t];
wx = [wx x];
end
x(t)
Metoda Eulera. Algorytm.
input t0, tf , x0, dt, tout
function [t, x] Ź Integrator (t, x, h , tend)
t Ź t0
while t < tend do
x Ź x0
if tend t < h then
m Ź 1
h Ź tend t
tpm Ź t
end
xpm Ź x
[t, x] Ź Euler (t, x, h)
while t < tf do
end
tend Ź t + tout
if tend > tf then
tend Ź tf
function [t, x] Ź Euler (t, x, h)
end
dxdt Ź Pochodna (t, x)
h Ź dt
x Ź x + dxdt * h
[t, x] Ź Integrator (t, x, h, tend)
t Ź t + h
m Ź m + 1
tpm Ź t
function dxdt Ź Pochodna (t, x)
dt krok
xpm Ź x
dxdt Ź ...
tout odstęp wyjścia
end do
tpm tablica czasu
output tpm, xpm
Nie wyprowadza
xpm tablica funkcji
wszystkich wyników
m licznik iteracji
Metoda Heuna
dx
= f (t, x)
dt
x0(t + h) = x(t) + f (t, x)h
nachylenie1 = f (t, x)
nachylenie2 = f (t + h, x0)
f (t, x) + f (t + h, x0)
nachylenieśr =
2
x(t + h) x(t) + f (t, x)h
metoda Eulera
Predyktor
Błąd lokalny metody Heuna
x0(t + h) = x(t) + f (t, x)h
óó
f (t +qh)
Korektor
Et = - h3 = O(h3)
f (t, x) + f (t + h, x0(t + h))
12
x(t + h) x(t) + h
Błąd globalny metody O(h2)
2
Metoda Heuna. Przykład. Algorytm
dx
8
= -2t3 +12t2 - 20t + 8.5
analityczne
dt
7
Heun h=0.5
Euler h=0.5
6
W tym przykładzie
x(0) =1
nie liczymy x0
5
(x nie występuje po
t = 0;
4
prawej stronie r.r.).
tf = 4;
3
h = input ('h = ');
x = 1;
2
wt = t;
1
wx = x; 0 1 2 3 4
t
while t < tf
dxdt1 = -2*t^3 + 12*t^2 - 20*t + 8.5;
function [t, x] Ź Heun (t, x, h)
t = t + h;
dxdt1 Ź Pochodna (t, x)
dxdt2 = -2*t^3 + 12*t^2 - 20*t + 8.5;
x0 Ź x + dxdt1 * h
x = x + (dxdt1 + dxdt2)/2*h;
dxdt2 Ź Pochodna (t + h, x0)
wt = [wt t];
nachylenie Ź (dxdt1 + dxdt2) / 2
Jeden krok
wx = [wx x];
x Ź x + nachylenie * h
metody
end
t Ź t + h
x(t)
Metoda Heuna. Iteracja korektora
Predyktor
function [t, x] Ź HeunIter (t, x, h)
x0(t + h) = x(t) + f (t, x)h
es Ź 0.0001
maxit Ź 20
Korektor
f (t, x) + f (t + h, x0(t + h))
dxdt1 Ź Pochodna (t, x)
x(t + h) x(t) + h
x0 Ź x + dxdt1 * h
2
for i Ź 1 to maxit do
x(t + h) występuje po obu stronach równania,
dxdt2 Ź Pochodna (t + h, x0)
więc
xe Ź x + (dxdt1 + dxdt2) / 2 * h
f (t, x) + f (t + h, xi-1(t + h))
if | xe - x0 | Ł es * |xe| then
xi (t + h) x(t) + h
break
2
end
Kryterium zakończenia iteracji
x0 Ź xe
end
xi (t + h) - xi-1(t + h)
x Ź xe
ea =
t Ź t + h
xi (t + h)
Jeden krok
metody
Metoda Heuna. Przykład 2
100
dx
= 4e0.8t - 0.5x
analitycze
dt
1 iteracja. korektora
80
15 iteracji korektora
x(0) = 2
60
Rozwiązanie analityczne:
40
4
x(t) = (e0.8t - e-0.5t ) + 2e-0.5t
1.3
20
Liczba iteracji
korektora
0
0 1 2 3 4
t
t x x błąd x błąd
(dokł.) 1 iter. wzg. % 15 iter. wzg. %
0 2.0000 2.0000 0 2.0000 0
1.0000 6.1946 6.7011 8.1756 6.3609 2.6835
2.0000 14.8439 16.3198 9.9425 15.3022 3.0876
3.0000 33.6772 37.1992 10.4584 34.7433 3.1657
4.0000 75.3390 83.3378 10.6171 77.7351 3.1805
x(t)
Metoda punktu środkowego
dx
= f (t, x)
dt
h
h
x(t + ) x(t) + f (t, x)
2
2
hh
nachylenie = f (t + , x(t + ))
22
Błąd lokalny metody O(h3)
h h
x(t + h) x(t) + f (t + , x(t + ))h
Błąd globalny metody O(h2)
2 2
Metoda punktu środkowego. Przykład. Algorytm
function [t, x] Ź Midpoint (t, x, h)
dx
dxdt Ź Pochodna (t, x)
= -2t3 +12t2 - 20t + 8.5
xh2 Ź x + dxdt * h / 2
dt
dxdth2 Ź Pochodna (t + h/2, xh2)
x(0) =1
Jeden krok x Ź x + dxdth2 * h
metody t Ź t + h
8
analityczne
7
punktu srodkowego
Eulera
6
5
4
3
2
h = 0.5
1
0 1 2 3 4
t
x(t)
Metoda Rungego-Kutty drugiego rzędu (1)
Metody Rungego-Kutty osiągają dokładność wzoru Taylora bez konieczności
liczenia pochodnych wyższego rzędu.
Metoda drugiego rzędu jest dana zależnością
dx
= f (t, x)
x(t + h) = x(t) + (a1k1 + a2k2)h
dt
gdzie
nachylenie
k1 = f (t, x)
k2 = f (t + ph, x + qk1h)
Porównując wzór metody z rozwinięciem x(t+h) wg wzoru Taylora dostajemy
a1 + a2 =1
1
3 równania
a2 p =
2
4 zmienne
1
Trzeba zatem założyć wartość
a2q =
2
jednej z niewiadomych.
Metoda Rungego-Kutty drugiego rzędu (2)
Niech a2 = 1/2. Wtedy a1 = 1/2, zaś p = q = 1. Zatem
1
x(t + h) = x(t) + (1 k1 + k2)h
2 2
Jest to metoda Heuna
gdzie
nachylenie
z pojedynczym korektorem
k1 = f (t, x)
k2 = f (t + h, x + k1h)
Niech a2 = 1. Wtedy a1 = 0, zaś p = q = 1/2. Zatem
x(t + h) = x(t) + k2h
Jest to metoda punktu
gdzie
nachylenie
środkowego
k1 = f (t, x)
1 1
k2 = f (t + h, x + k1h)
2 2
Porównanie metod
8
analityczne
h = 0.5
punktu srodk.
7
Eulera
Heuna
6
5
4
3
2
Metody
Metoda
1
0 1 2 3 4
2-go rzędu
1-go rzędu
t
Euler Heun Punkt
środkowy
Błąd lokalny metody O(h2) O(h3) O(h3)
Błąd globalny metody O(h) O(h2) O(h2)
x(t)
Metoda Rungego-Kutty trzeciego rzędu
Metoda trzeciego rzędu jest dana zależnością
1
dx
x(t + h) = x(t) + (k1 + 4k2 + k3)h
= f (t, x)
6
dt
gdzie
nachylenie
k1 = f (t, x)
1 1
k2 = f (t + h, x + k1h)
2 2
1
k3 = f (t + h, x - k1h +2k2h)
2
Błąd lokalny metody O(h4)
Metoda Rungego-Kutty czwartego rzędu
Metoda czwartego rzędu jest dana zależnością
dx
= f (t, x)
dt
1
x(t + h) = x(t) + (k1 + 2k2 + 2k3 + k4)h
6
gdzie
nachylenie
k1 = f (t, x)
1 1
k2 = f (t + h, x + k1h)
2 2
1 1
k3 = f (t + h, x + k2h)
2 2
k4 = f (t + h, x +k3h)
Błąd lokalny metody O(h5)
Metoda Rungego-Kutty czwartego rzędu. Algorytm
dx
= f (t, x)
function [t, x] Ź RK4 (t, x, h)
dt
k1 = f (t, x)
k1 Ź f (t, x)
xm Ź x + k1 * h / 2
1 1
k2 = f (t + h, x + k1h)
2 2
k2 Ź f (t + h/2, xm)
1 1
k3 = f (t + h, x + k2h)
xm Ź x + k2 * h / 2
2 2
k3 Ź f (t + h/2, xm)
k4 = f (t + h, x +k3h)
xe Ź x + k3 * h
k4 Ź f (t + h, xe)
nachylenie = (k1+2k2+2k3+k4)/6
1
x Ź x + nachylenie * h
x(t + h) = x(t) + (k1 + 2k2 + 2k3 + k4)h
6
t Ź t + h
Jeden krok
metody
Metoda Rungego-Kutty czwartego rzędu. Przykłady
dx
dx
= 4e0.8t - 0.5x
= -2t3 +12t2 - 20t + 8.5
dt
dt
x(0) = 2
x(0) =1
5
80
analityczne
4.5
70
Runge-Kutta 4
4
60
3.5
50
3
40
2.5
30
2
20
analityczne
1.5
10
Runge-Kutta 4
1
0
0 1 2 3 4
0 1 2 3 4
h = 0.5 h = 1
Problemy z rozwiązaniami, które wykazują gwałtowną zmianę
2
dx
+ 0.6x =10e-(t-2) /(20.0752 )
dt
x(0) = 0.5
2
1.8
h = 0.1
1.6
1.4
Zmniejszenie kroku
1.2
1
0.8
Zwiększenie liczby
obliczanych wartości
0.6
funkcji
0.4
0.2
Alternatywa: krok
0
o zmiennej długości
0 0.5 1 1.5 2 2.5 3 3.5 4
Metoda adaptacyjna Rungego-Kutty-Fehlberga (1)
Metody włączenia zmiany długości kroku:
1. Lokalny błąd jest szacowany w oparciu o tą samą metodę RK użytą
dla dwóch kroków o różnych długościach
2. Lokalny błąd obcięcia jest szacowany jako różnica pomiędzy dwoma
przewidywaniami opartymi o dwie metody RK o różnych rzędach
Liczba wartości funkcji obl. w jednym kroku 1 2 3 4 5 6 7 8
Maksymalny rząd metody Rungego-Kutty 1 2 3 4 4 5 6 6
Fehlberg
wyższy rząd metody w każdym kroku więcej wartości funkcji f do obliczeń
Tak dobrał wartości parametry metody RK, aby w obu metodach występowały
te same wartości. Dzięki temu wystarczy 6 wartości zamiast 10. Metoda RKF
jest 5-tego rzędu.
Metoda adaptacyjna Rungego-Kutty-Fehlberga (2)
Przybliżenie wartości x(t+h) odpowiednio czwartego i piątego rzędu:
37 250 125 512
ć
x(t + h) = x(t) + k1 + k3 + k4 + k5 h
378 621 594 1771
Ł ł
2825 18575 13525 277 1
ć
x(t + h) = x(t) + k1 + k3 + k4 + k5 + k6 h
27648 48384 55296 14336 4
Ł ł
gdzie
k1 = f (t, x)
Różnica daje oszacowanie
1 1
k2 = f (t + h, x + k1h)
5 5
aktualnego błędu
3 3 9
k3 = f (t + h, x + k1h + k2h)
10 40 40
dx
3 3 9 6
= f (t, x)
k4 = f (t + h, x +10 k1h -10 k2h + k3h)
5 5
dt
5 70 35
k5 = f (t + h, x -11 k1h + k2h - k3h + k4h)
54 2 27 27
7 1631 575 44275 253
k6 = f (t + h, x + k1h +175 k2h +13824k3h + k4h + k5h)
8 55296 512 110592 4096
Metoda adaptacyjna Rungego-Kutty-Fehlberga (3)
Zmiana wielkości kroku:
Sposób I Sposób II
a
hnowy = 2hobecny
Dwymagany
hnowy = haktualny
Daktualny
gdy Daktualny Ł aDwymagany
gdzie
gdzie 0 < a <1
Dwymagany żądana dokładność
Daktualny błąd aktualny
haktualny
hnowy =
a parametr
2
gdy Daktualny > Dwymagany
a = 0.2 gdy krok zwiększamy
tj. gdy Daktualny Ł Dwymagany
a = 0.25 gdy krok zmniejszamy
tj. gdy Daktualny > Dwymagany
Metoda adaptacyjna Rungego-Kutty-Fehlberga (4)
2
2
1.8
dx
+ 0.6x =10e-(t-2) /(20.0752 )
RKF
1.6
dt
1.4
x(0) = 0.5
1.2
1
0.8
RK4
0.6
2
0.4
1.8
RK4
0.2
1.6
0
1.4
0 0.5 1 1.5 2 2.5 3 3.5 4
1.2
RKF
1
0.8
0.6
0.4
0.2
0
0 0.5 1 1.5 2 2.5 3 3.5 4
Wyszukiwarka
Podobne podstrony:
Wyklad studport 8
wyklad 7 zap i, 11 2013
socjo wykład z 26 11
5 Analiza systemowa wykłady PDF 11 z numeracją
wyklad 8 zap i, 11 2013
Techniki negocjacji i mediacji w administracji wykłady 05 11 2013
Analiza Wykład 8 (25 11 10)
KPC Wykład (7) 13 11 2012
Wyklad studport 9
KPC Wykład (6) 06 11 2012
Enzymologia wyklad 30 11
Materiały do wykładu 7 (18 11 2011)
Wyklad WektoryMacierze 11 08
Wyklad AnalizaMat 11 08
więcej podobnych podstron