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 Pa��s
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 4m�L dh D�p
+� +� h =�
dh
2g dt 2r�g
dt2 r�g R2
h(0) =� 0, =� 0
dt
t=�0
Manometr (2)
2
L d h 4m�L dh D�p
+� +� h =�
dx1
2g dt 2r�g
dt2 r�g R2
=� x2
dt
2
d h 8m� dh 2g 1
dx2
+� +� h =� D�p
=� cD�p -� ax2 -� bx1
dt L r�L
dt2 r� R2
dt
8m� 2g 1
Warunki początkowe
a =� , b =� , c =�
L r�L
r� R2
x1(0) =� 0, x2(0) =� 0
2
d h dh
+� a +� bh =� cD�p
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 +�q�h)
ó�ó� 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 +�q�h)
Rn =� hn+�1 gdzie 0 <�q� <�1
(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 +�q�h)
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 +�q�h)
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
e�a =�
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) /(2��0.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
D�wymagany
hnowy =� haktualny
D�aktualny
gdy D�aktualny Ł� a��D�wymagany
gdzie
gdzie 0 < a <1
D�wymagany  żądana dokładność
D�aktualny  błąd aktualny
haktualny
hnowy =�
a�  parametr
2
gdy D�aktualny > D�wymagany
a� = 0.2 gdy krok zwiększamy
tj. gdy D�aktualny Ł� D�wymagany
a� = 0.25 gdy krok zmniejszamy
tj. gdy D�aktualny > D�wymagany
Metoda adaptacyjna Rungego-Kutty-Fehlberga (4)
2
2
1.8
dx
+� 0.6x =�10e-�(t-�2) /(2��0.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