Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Równania różniczkowe zwyczajne
Zagadnienie początkowe
dy
= f ( y(x), x), y(a) = c
dx
Szukamy różniczkowalnej y(x) spełniającej powyższe warunki.
Warunek Lipschitza:
"x "[a,b], y1, y2 " R "L > 0 f ( y1, x) - f ( y2, x) d" L y1 - y2
f : [a,b]� R R
Jeżeli ciągła funkcja spełnia warunek Lipschitza:
to zagadnienie początkowe ma dokładnie jedno różniczkowalne w
sposób ciągły rozwiązanie y:[a,b]R.
W9-1
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
rozwiązanie numeryczne: na przedziale [a,b]
punkty: a=x0, x1, x2, ....., xn=b, xi-xi-1=hi
SCHEMAT
SCHEMAT
RÓŻNICOWY
RÓŻNICOWY
wartości przybliżone: y0, y1, y2, ....., yn
ANALIZA BA

DU
?
wartości dokładne: y(x0), y(x1), .......,y(xn)
Błędy
W9-2
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Wybrane schematy różnicowe jednokrokowe
Metoda EULERA
yn+1 - yn yn - yn+1
dy
= f ( y( x ), x ), H" f ( yn , xn ) = f ( yn+1 , xn+1 )
dx h - h
yn+1 = yn + f ( yn , xn )h yn+1 = yn + f ( yn+1 , xn+1 )h
otwarta (jawna) zamknięta (niejawna)
Zmodyfikowana metoda Eulera
pół kroku w kierunku pochodnej - poprawienie pochodnej - cały krok w
poprawionym kierunku
h h
fn = f ( yn , xn ) y = yn + fn f = f ( y , xn + )
1 1 1
n+ n+ n+
2 2
2 2 2
yn+1 = yn + f h
1
n+
2
Schemat jednokrokowy
yn+1 = yn + Ś ( yn+1 , yn , xn ,h )h - niejawny (zamknięty)
f
yn+1 = yn + Ś ( yn , xn ,h )h - jawny (otwarty)
f
W9-3
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Schematy Rungego-Kutty r-poziomowy (etapowy) schemat R-K
r
Ś ( yn , xn ,h ) =
f "c Ki( yn , xn ,h )
i
i =1
r r
Ki ( y, x,h ) = f ( y + h i = 1,2,...,r
"b K , x + h"b ),
ij j ij
j =1 j =1
otwarty (jawny): bij = 0 j e" i
i -1 i -1
K1 = f ( y, x ), Ki ( y, x,h ) = f ( y + h i = 2,...,r
"b K , x + h"b ),
ij j ij
j =1 j =1
K1 = f ( yn , xn ),
1 1
K2 = f ( yn + K1 , xn + h ),
2 2
RK4:
1 1
K3 = f ( yn + K2 , xn + h ),
2 2
K4 = f ( yn + K3 , xn + h ),
1
yn+1 = yn + h( K1 + 2K2 + 2K3 + K4 )
6
W9-4
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Zbieżność metod jednokrokowych i sterowanie długością kroku
yn+1 = yn + Ś ( yn , xn ,h )h
f
yn+1 = y( xn ) + Ś ( y( xn ), xn ,h )h
f
rn+1 = y( xn + h ) - [ y( xn ) + Ś ( y( xn ),xn ,h )h] - błąd lokalny
f
Metoda rządu p:
p+1
p+1 p+2
rn+1( h ) = �( y( xn ), xn )h + O( h ) (rn+1(h) < Ch )
metoda p
Eulera 1
zmodyfikowana Eulera 2
RK2,3,4 2,3,4
RKm m=5,6,.. pJeżeli metoda jednokrokowa jest rzędu p i spełnia warunek Lipschitza
~
Ś (y, x,h) - Ś (~, x,h) d" LŚ y - y
y
to błąd globalny można oszacować przez
f f
p
().
max yl - y(xl ) d" Chmax , hmax = max xl+1 - xl
l=0,1, ,n l=0,1, ,n-1
W9-5
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
xn , y( xn ) yn+1 xn , y( xn ) y
1 1
h h
h n+ +
2 2
2 2
p+1
h
p+1
y( xn + h ) - yn+1 = �h + ... y( xn + h ) - y = 2��# ś# + ...
ś# ź#
1 1
n+ +
2
�# #
2 2
p
2
p+1
ERR=�h H" ( y - yn+1 )
1 1
p
n+ +
2 - 1
2 2
16
dla p=4 : ERR= ( y - yn+1 )
1 1
n+ +
15
2 2
16 y - yn+1
1 1
n+ +
2 2
ERR < ymaxRELREER+ABSERR �! yn+1 :=
15
ERR > ymaxRELREER+ABSERR �! zmniejszyć h
W9-6
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
inaczej
yn + yn+1
ERR < RELREER+ABSERR lub
2
ERR
ETOL= < 1
yn + yn+1
RELERR + ABSERR
2
5
szukamy ą, by ETOL(ąh ) H" 1 : ETOL(ąh ) H" ą ETOL( h )
1
ą = 0.9, hnowe = ąh
5
ETOL( h )
�# ś#
0.9
ś# ź#
Jeżeli ETOL( h ) < 1, to ą = minś# 5 ,5ź#
ETOL( h )
�# #
Metoda wymaga dla RK4 4+7=11 obliczeń prawej strony.
W9-7
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Algorytm Rungego-Kutty-Fehlberga stosuje dwa schematy
Rungego-Kutty m+1 i m etapowy z odpowiednio dobranymi
współczynnikami. Schemat m-etapowy jest rzędu p a schemat m+1
etapowy jest rzędu p+1, a współczynniki Ki są jednakowe.
m +1
yn+1 = yn + Ś ( yn , xn ,h )h = yn +
f "c Ki ( yn , xn ,h )
i
i =1
m
~n+1 = yn + Ś f ( yn , xn ,h )h = yn + ~iKi ( yn , xn ,h ) ~
y cm+1 = 0
"c
i =1
m +1
~
Wtedy ERR( h ) = h
"( ci - ci )
i =1
W9-8
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Metody wielokrokowe
Metody Adamsa: Adamsa-Bashfortha (jawne), Adamsa  Moultona
(niejawne)
xn+1
y( xn+1 ) = y( xn ) + f ( y( x ), x )dx
+"
xn
przybliżamy funkcję podcałkową wielomianem interpolacyjnym
Lagrange a w węzłach
i = -1,0,1,...k metoda niejawna
ż#
( xn-i , f ( yn-i , xn-i ))
�#
i = 0,1,...k metoda jawna
�#
całkujemy wielomian i otrzymujemy wzory postaci
k k
~
~j
yn+1 = + h � f ( yn- j , xn- j ) dla metody jawnej
"ą yn- j " j
j =0 j =0
k k
yn+1 = yn- j + h � f ( yn- j , xn- j ) + h�-1 f ( yn+1 , xn+1 ) dla
"ą j " j
j =0 j =0
metody niejawnej
W9-9
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Metoda wstecznego różniczkowania
przybliżamy rozwiązanie y(x) wielomianem interpolacyjnym W(x)
zbudowanym na węzłach
i = -1,0,1,...k metoda niejawna
ż#
( xn-i , yn-i )
�#
i = 0,1,...k metoda jawna
�#
obliczamy pochodną tego wielomianu W (x) , która przybliża
pochodną rozwiązania,
z równości
2
y (xn ) = f (yn, xn) H" W '(xn) otrzymujemy wzór
k
~
~j
yn+1 = + h�0 f (yn, xn) dla metody jawnej
"ą yn- j
j=0
2
a z y (xn+1) = f (yn+1, xn+1) H" W '(xn+1)
k
yn+1 = yn- + h�-1 f (yn+1, xn+1) dla metody niejawnej
"ą j j
j =0
W9-10
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Stosowanie metod niejawnych
k k
yn+1 = yn- j + h � f ( yn- j , xn- j ) + h�-1 f ( yn+1 , xn+1 )
"ą j " j
j =0 j =0
jest nieliniowym równaniem algebraicznym względem yn+1 i jest
rozwiązywane metodą iteracyjną (najczęściej iteracji prostej):
k k
[ i [ i
yn+] = yn- j + h , xn- j ) + h�-1 f ( yn+-1 ] , xn+1 ), i=1,2,...
1 "ą j "� f ( yn- j 1
j
j =0 j =0
Wymaga ona punktu startowego, np. z metody Eulera
[ 0 ]
yn+1 = yn + hf ( yn , xn ). Uzyskanie dokładnego rozwiązania wymaga
wtedy wielu iteracji.
Metody typu PREDYKTOR-KOREKTOR wyznaczają przybliżenie
początkowe jawną metodą wielokrokową o takiej samej liczbie
kroków, a następnie stosują kilka iteracji rozwiązujących równanie
nieliniowe.
Algorytm GEAR A
W9-11
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Liniowe układy równań różniczkowych
d
y( x ) = Ay( x ), y( 0 ) = y0
dx
gdzie A jest macierzą nxn o różnych wartościach własnych
i i = 1,2,...,n

vi i = 1,2,...,n
i wektorach własnych

Avi = ivi , vi `" 0
Wartość własna i wektor własny:
( - A vi = 0 det Ii - A = 0
) ( )
Ii
czyli wartości własne są pierwiastkami równania
( )
det Is - A = 0
wielomian stopnia n (wielomian charakterystyczny A)
W9-12
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Niech:

[
V = v1 v2 vn ], Vz( x ) = y( x )
Wtedy:
1 0 0
Ą# ń#
ó# Ą#
0 2 0
ó# Ą#
� =
ó# Ą#
AV = V�
ó#
0 0 nĄ#
Ł# Ś#
-1
V AV = �
W9-13
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Niech:

Vz( x ) = y( x )
d
-1 -1
V z( x ) = AVz( x ), z( 0 ) = V y( 0 ) = V y0 = z0
dx
d
-1 -1 -1
z( x ) = V AVz( x ), z( 0 ) = V y( 0 ) = V y0 = z0
dx
1 0 0
Ą# ń#
ó# Ą#
0 2 0
d
ó# Ą#z( x ), z( 0 ) = V -1y( 0 ) = V -1y0 = z0
z( x ) =
ó# Ą#

dx
ó#
0 0 nĄ#
Ł# Ś#
W9-14
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
d
zi( x ) = izi( x ), zi( 0 ) = zi0 i = 1,2,...,n
dx
i
zi( x ) = e xzi0 i = 1,2,...,n
1
Ą# ń#
e x 0 0
ó# Ą#
2
0 e x 0

ó# Ą#z( 0 ),
z( x ) =
ó# Ą#

ó# Ą#
n
0 0 e x Ś#
ó# Ą#
Ł#
W9-15
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
1
Ą# ń#
e x 0 0
ó# Ą#
2
0 e x 0

ó# Ą#z( 0 ),
Vz( x ) = V
ó# Ą#

ó# Ą#
n
0 0 e x Ś#
ó# Ą#
Ł#
T ń#
Ą#
w1
1
Ą# ń#
e x 0 0
T
ó# Ą#
ó# Ą#
2
-1
2
ó#w Ą#
0 e x 0
V :=
ó# Ą#V -1y( 0 ),
y( x ) = V
ó# Ą#

ó# Ą#

ó# Ą#
T
ó# Ą#
n
0 0 e x Ś# ó#
ó# Ą#
n
Ł#w Ą#
Ś#
Ł#
W9-16
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
n

i x
y( x ) =
"e viwiT y( 0 ),
i=1
W9-17
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Stabilność i sztywność
dy
Rozważmy układ n równań = Ay( x ), y( 0 ) = y0
dx
n
x
Ax
i
rozwiązanie: y( x ) = e y0 =
"c e dąży do 0 dla Re i < 0
i
i =1
Metodą Eulera
yn+1 = yn + hAyn = ( I + hA )yn: y1 = ( I + hA )y0,
y2 = ( I + hA )y1 = ( I + hA )2 y0
y3 = ( I + hA )y2 = ( I + hA )3 y0
.............................................
lim yn = 0 �! 1 + ih < 1 i = 1,...,n
n"
Obszar stabilności schematu różnicowego: zbiór wartości
zespolonych ih, dla których wszystkie rozwiązania zadania
testowego są ograniczone dla n ". Jeżeli obszar stabilności
zawiera punkt 0, to metodę nazywamy stabilną.
Metoda A-stabilna A(a)-stabilna , A(ą) stabilna
W9-18
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
W9-19
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
W9-20
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
W9-21
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
W9-22
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
W9-23
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
W9-24
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
W9-25
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Obszary stabilności absolutnej
Adamsa-Bashfortha
l. kroków 1 2 3 4
rząd 1 2 3 4
q* -2 -1 -6/11 -3/10
Adamsa-Moultona
l. kroków 1 2 3 4
rząd 2 3 4 5
q* -6 -3 -90/49
-"
W9-26
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Wstecznego różniczkowania
k 1 2 3 4 5 6
a 0 0 -0.1 -0.7 -2.4 -6.1
90 90 88 73 51 18
ą
W9-27
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
W9-28
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Rungego-Kutty
p=m 1 2 3 4
Lewa -2 -2 -2.51 -2.78
granica
obszaru
stabilności
Żadna metoda jawna nie jest A-
stabilna A(a)-stabilna lub A(ą)-
stabilna dla żadnych a i ą.
Każda metoda jawna ma ograniczony obszar stabilności.
W9-29
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
W9-30
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
ode45 is based on an explicit Runge-Kutta (4,5) formula, the
Dormand-Prince pair. In general, ode45 is the best solver to apply
as a first try for most problems. For this reason, ode45 is the default
solver used by Simulink for models with continuous states.
W9-31
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
ode23 is also based on an explicit Runge-Kutta (2,3) pair of Bogacki
and Shampine. It can be more efficient than ode45 at crude
tolerances and in the presence of mild stiffness.
ode113 is a variable-order Adams-Bashforth-Moulton PECE solver.
It can be more efficient than ode45 at stringent tolerances.
ode15s is a variable order solver based on the numerical
differentiation formulas (NDFs). These are related to but are more
efficient than the backward differentiation formulas, BDFs. If you
suspect that a problem is stiff, or if ode45 failed or was very
inefficient, try ode15s.
ode23s is based on a modified Rosenbrock formula of order 2.
Because it is a one-step solver, it can be more efficient than ode15s
at crude tolerances. It can solve some kinds of stiff problems for
which ode15s is not effective.
ode23t is an implementation of the trapezoidal rule using a "free"
interpolant. Use this solver if the problem is only moderately stiff
and you need a solution without numerical damping.
W9-32
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
ode23tb is an implementation of TR-BDF2, an implicit Runge-Kutta
formula with a first stage that is a trapezoidal rule step and a
second stage that is a backward differentiation formula of order
two. By construction, the same iteration matrix is used in evaluating
both stages. Like ode23s, this solver can be more efficient than
ode15s at crude tolerances.
Note: For a stiff problem, solutions can change on a time scale that is
very short compared to the interval of integration, but the solution of
interest changes on a much longer time scale. Methods not designed for
stiff problems are ineffective on intervals where the solution changes
slowly because they use time steps small enough to resolve the fastest
possible change. Jacobian matrices are generated numerically for ode15s
and ode23s. For more information, see Shampine, L. F., Numerical
Solution of Ordinary Differential Equations, Chapman & Hall, 1994.
W9-33
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Przykłady
Rozwiązanie wolnozmienne
df (t) f (t) 3
ś#
= 0.1 f (t)�#1- f (0) = 1, 0 < t < 50, f (t) =
ś# ź#
-t
dt 3 10
�# #
1+ 2e
Błąd w funkcji RelTol
-2
10
RKF23
RKF45
-4
Adams PC
10
Gear WR
-6
10
-8
10
-10
10
-12
10
-14
10
-12 -10 -8 -6 -4 -2
10 10 10 10 10 10
W9-34
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
Kroki w funkcji RelTol
4
10
RKF23
RKF45
Adams PC
Gear WR
3
10
2
10
1
10
-12 -10 -8 -6 -4 -2
10 10 10 10 10 10
W9-35
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne w Inżynierii wykład 9
LW F w funkcji RelTol
5
10
RKF23
RKF45
Adams PC
Gear WR
4
10
3
10
2
10
1
10
-12 -10 -8 -6 -4 -2
10 10 10 10 10 10
W9-36
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne ET3 wykład 6
Równanie sztywne
dy (t) df (t)
= -a(y(t) - f (t))+ y(0) = 2, a = 500
dt dt
3
- at
0 < t < 50 , f (t) = , y(t) = e + f (t)
- t
10
1 + 2e
3
2.8
2.6
2.4
2.2
2
1.8
1.6
1.4
1.2
1
0 5 10 15 20 25 30 35 40 45 50
W6-37
Instytut Automatyki Politechniki A
ódzkiej - Metody Numeryczne ET3 wykład 6
LW F w funkcji RelTol
5
10
RKF23
RKF45
Adams PC
Gear W R
4
10
3
10
2
10
1
10
-6 -5 -4 -3 -2
10 10 10 10 10
W6-38