Rozwiązywanie równań różniczkowych 1
11 Rozwiązywanie równań różniczkowych
Przykład 1
Należy obliczyć przebieg czasowy u(t), t " [t0, t0 + T ] w układzie elektrycznym (ozn.RC):
R
Równanie dynamiki tego układu:
dq(t) e(t) - u(t)
=
dt R
e(t) q(t) u(t)
q(t0) = Q(u(t0))
Dla staÅ‚ej, liniowej pojemnoÅ›ci Q(u) = C · u, czyli
Å„Å‚q(t) = Q(u(t)).
òÅ‚
0 dla t < t0
JeÅ›li e(t) = E0 · 1(t - t0), gdzie 1(t - t0) = 0, to
ół
1 dla t t0
0
istnieje rozwiÄ…zanie analityczne: u(t) = E0 · (1 - e-(t-t )/Ä)1(t - t0).
W ogólnym przypadku można jedynie uzyskać numeryczne przybliżenie rozwiązania.
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 2
Standardowa postać układu równań różniczkowych
dyi(t)
= fi(t, y(t)), i = 1, . . . , m, gdzie
dt
y(t) = [y1(t); y2(t); . . . ym(t)]
Zapis zwarty układu równań różniczkowych:
y (t) = f(t, y(t)), gdzie
f(t, y) = [f1(t, y); f2(t, y); . . . ; fm(t, y)]
Uwaga: Równanie różniczkowe rzędu n > 1 można sprowadzić do układu n równań
różniczkowych pierwszego rzędu.
Przykład 2
Równanie oscylatora harmonicznego: x (t) = -x(t) można sprowadzić do układu:
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
y1(t) y2(t)
ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚, przy czym x(t) a" y1(t)
=
y2(t) -y1(t)
y (t) f(t,y(t))
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 3
11.1 Sformułowanie zadania
Ogólne sformułowania problemu rozwiązywania układu równań różniczkowych
" Zagadnienie początkowe (Cauchy ego): obliczyć przebieg y(t) w przedziale
czasu t " [t0, tF ] (gdzie tF = t0+T ) znając wektor rozwiązań początkowych y(t0) = y0
y (t) = f(t, y), dla t " [t0, tF ], gdzie y(t0) = y0
" Zagadnienia brzegowe: obliczyć przebieg y(t) w przedziale czasu t " [t0, tF ] spełniający
zadany warunek brzegowy postaci Ay(t0) + By(tF ) = C.
Np. dla zagadnienia prostego strzału trzeba znalez
ć takie składowe wektora y(t0)
(parametry  strzału ), dla których rozwiązanie przechodzi w pewnej chwili tF przez zadany
punkt yF ( cel strzału )
y (t) = f(t, y), dla t " [t0, tF ], gdzie y(tF ) = yF
Dalej zajmujemy siÄ™ jedynie metodami dyskretnymi rozwiÄ…zywania zagadnienia poczÄ…tkowego,
tzn. przybliżenie rozwiązania będzie wyznaczane w skończonej liczbie punktów
t0 < t1 < t2 < · · · < tN = tF .
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 4
Własności zadania początkowego
Jeżeli funkcje fi(t, y), i = 1, . . . , n
" są ciągłe na zbiorze S = {t, y) : t0 t tF , yi " R} i
" spełniają na tym zbiorze warunek Lipschitza ze względu na y, tj.
"L > 0: "t " [t0, tF ] i dowolnych y1, y2 " Rn:
n
|fk(t, y1) - fk(t, y2)| L y1 - y2 , k = 1, . . . , n
i=1
to zagadnienie początkowe przy dowolnym y(t0) = y0 " R ma dokładnie jedno rozwiązanie
1
y(t) klasy C[t ,tF ].
0
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 5
Zbieżność metod dyskretnych rozwiązywania zagadnienia początkowego
Jeżeli punkty ti, i = 0, . . . N sÄ… równoodlegÅ‚e, tzn. ti = t0 + i · h (gdzie
h = T/N, T = tF - t0 jest krokiem całkowania), to:
" Błąd przybliżenia yi rozwiązania y(t) w punkcie ti jest zdefiniowany jako:
ei = yi - y(ti)
" Metoda rozwiÄ…zania jest zbieżna, jeÅ›li "µ > 0, "h0 = h0(µ) takie, że "h " (0, h0)
zachodzi: a" maxi=0,...N(h) ei < µ.
g
jest globalnym błędem metody dla danego zadania i długości kroku h.
g
" Największą liczbę całkowitą p taką, że dla dostatecznie małych h zachodzi:
a" max ei < c · hp
g
i=0,...N(h)
gdzie c jest stałą niezależną od h, nazywamy rzędem metody rozwiązywania zagadnienia
poczÄ…tkowego.
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 6
11.2 Rozwiązywanie równań różniczkowych w programie MATLAB
Przybliżone rozwiązanie układu równań różniczkowych:
y (t) = f(t, y), y(t0) = y0
w przedziale czasu [t0, tf] można uzyskać za pomocą wywołania o postaci
[ T , Y ] = s o l v e r ( odefun , [ t0 , t f ] , y0 , o pt io ns )
gdziesolverjest nazwą metody całkowania:ode45, ode23, ode113,
ode15s, ode23s, ode23t, ode23tb(patrz np.doc ode45).
MacierzYzawiera przybliżenia wektora rozwiązań y w chwilach czasowych zawartych w
wektorzeT.
M-funkcjaodefun, która wyznacza wartości funkcji f(t, y) dla zadanych wartości
argumentów t i y, ma postać
function dydt = odefun ( t , y )
n=length ( y ) ;
dydt ( 1 ) = . . . ; dydt ( 2 ) = . . . ; dydt ( n)= . . . . ;
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 7
Przykład 3
R R
1 2
Modelem matematycznym układu elektrycznego
(ozn. R2C2) pokazanego na rysunku Ò!
e(t)
u1 u2
jest następujący układ równań różniczkowych:
C C
1 2
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
e(t)
R1+R 1
u 1(t) -R1R2C2 R2C1 ûÅ‚ ðÅ‚ u1(t)
R1C1
1
ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚ ðÅ‚ ûÅ‚
= +
1
u 2(t) -R21C2 u2(t) 0
R2C2
u(t) f(t,u(t))
Powyższy układ można uzyskać przez przekształcenie poniższych równań, wynikających z pierwszego i drugiego prawa
Kirchoffa oraz równań dynamiki kondensatorów:
ie = iC1 + iC2
e = u1 + ieR1, u1 = u2 + i2R2
iC1 = C1u 1, iC2 = C2u 2
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 8
Dla R1 = 1k&!, R2 = 1M&!, C1 = C2 = 1µ F, e(t) = 1(t)
Układ równań: Reprezentacja równań wMATLAB-ie:
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
function dudt=r2c2fun ( t , u )
3
-1001 1
ðÅ‚u1(t)ûÅ‚=ðÅ‚ ûÅ‚ ðÅ‚u1(t) ûÅ‚+ðÅ‚10 · 1(t) ûÅ‚
dudt (1,1)=-1001"u (1)+ u ( 2 ) ;
u 2(t) 1 -1 u2(t) 0
i f t >0 , dudt (1 ,1)= dudt (1 ,1)+1000; end
dudt (2 ,1)= u(1)-u ( 2 ) ;
A
RozwiÄ…zanie analityczne dla ui(0) = 0:
u1(t) = c1(1 - 1)e1t + c2(1 - 2)e2t + 1
u2(t) = c1e1t + c2e2t + 1, t 0
c1 = 2/(1 - 2), c2 = -1 - c1
" "
1 = -501 + 250001 H" -0.999, 2 = -501 - 250001 H" -1001
1, 2 są wartościami własnymi macierzy A.
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 9
opts=odeset ( RelTol ,1e-4, AbsTol ,1e-12);
RozwiÄ…zanie numeryczne uzyskane
[ t , y ]= ode23 r2c2fun , [ 0 , 0 . 1 ] , [ 0 ; 0 ] , opts ) ;
(
za pomocÄ… programuMATLAB
N=max( size ( t ) ) ;
y r e f =zeros (N, 2 ) ; y e r r =zeros (N, 2 ) ;
rozwiazania
1
c1=-0.5"(501/ sqrt (250001)+1); c2=-1-c1 ;
u1
lambda1=-501+sqrt (250001);
u2
0.5
lambda2=-501-sqrt (250001);
y r e f ( : , 2 ) = 1 + c1"exp ( lambda1" t )+ c2"exp ( lambda2" t ) ;
0
0 0.02 0.04 0.06 0.08 0.1
y r e f ( : , 1 ) = 1 + c1"(1+lambda1 )"exp ( lambda1" t ) + . . .
bledy wzgledne
-2
u1
10
c2"(1+lambda2 )"exp ( lambda2" t ) ;
u2
for i =2:N
-4
10
y e r r ( i , : ) = ( y ( i ,:) - y r e f ( i , : ) ) . / y r e f ( i , : ) ;
-6
end
10
0.02 0.04 0.06 0.08 0.1
subplot ( 3 1 1 ) ;
-3
krok calkowania
x 10
plot ( t , y ( : , 1 ) , - , t , y ( : , 2 ) , : ) ;
4
subplot ( 3 1 2 ) ;
2
semilogy ( t , abs ( y e r r ( : , 1 ) ) , - , t , abs ( y e r r ( : , 2 ) ) , : ) ;
subplot ( 3 1 3 ) ;
0
0 0.02 0.04 0.06 0.08 0.1
dt = d i f f ( t ) ;
t [s]
plot ( t ( 1 : length ( t )-1) , dt ) ;
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 10
11.3 Podstawy konstrukcji dyskretnych metod całkowania
" Metody jednokrokowe wykorzystują informację z jednego (ostatniego) kroku. Mają postać:
yn+1 = yn + hnÅš(tn, yn, hn), n = 0, 1, . . . , N - 1
" Liniowe metody wielokrokowe wykorzystują informację z K poprzednich kroków. Formuły
stałokrokowe mają postać:
K K
yn+1 = Ä…kyn+1-k + h ²k fn+1-k, n = 0, 1, . . . , N - 1
k=1 k=0
y n+1-k
gdzie fi = f(ti, y(ti)).
 JeÅ›li ²0 = 0, to formuÅ‚a jest jawna (otwarta, ekstrapolacyjna)
 JeÅ›li ² = 0, to formuÅ‚a jest uwikÅ‚ana (zamkniÄ™ta, interpolacyjna); nieznane yn+1

występuje po lewej i po prawej stronie formuły, więc wyznaczenie wartości tego wektora
wymaga rozwiÄ…zania ukÅ‚adu równaÅ„ algebraicznych (nieliniowych, jeÅ›li funkcja f(t, ·)
jest nieliniowa).
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 11
11.3.1 Metody jednokrokowe
Dla wyznaczenia przybliżonej wartości y(tn+1) metody jednokrokowe korzystają z wartości y i
y w jednej, poprzedniej chwili tn (i ew. w chwilach pośrednich).
" Jeśli wykorzystać przybliżenie różnicowe pochodnej z różnicą progresywną:
y(tn+1) - y(tn)
y (tn) H"
tn+1 - tn
to uzyskujemy metodę całkowania nazywaną metodą otwartą Eulera:
y(tn+1) = y(tn) + (tn+1 - tn)y (tn) = y(tn) +hn f(tn, y(tn))
hn yn
fn
Zapis zwięzły:
yn+1 = yn + hnfn, n = 0, 1, 2, ..., N - 1
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 12
" Jeśli wykorzystać przybliżenie różnicowe pochodnej z różnicą wsteczną:
y(tn+1) - y(tn)
y (tn+1) H"
tn+1 - tn
to uzyskujemy metodę całkowania nazywaną metodą zamkniętą Eulera:
y(tn+1) = y(tn) + (tn+1 - tn)y (tn+1) = y(tn) +hnf(tn+1, y(tn+1))
hn yn
Zapis zwięzły:
yn+1 = yn + hnfn+1, n = 0, 1, 2, ..., N - 1
Uwaga: niewiadoma n-tego kroku, t.j. yn+1, występuje po lewej i prawej stronie (ogólnie
nieliniowego) równania.
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 13
" Metody Rungego-Kuty mają postać:
R R
y(tn+1) = y(tn) + hn wrkr,n, gdzie kr,n = f(tn + hnar, yn + br,ÁkÁ,n)
r=1 Á=1
przy czym wr, ar, br,Á sÄ… parametrami (staÅ‚ymi) metody. Przypadki szczególne:
 Metoda Heuna:
hn
y(tn+1) = y(tn) + f(tn, yn) + f tn+1, y(tn) + hnf(tn, yn)
2
 Zmodyfikowana metoda Eulera:
hn hn
y(tn+1) = y(tn) + hnf tn + , yn + f(tn, y(tn))
2 2
  Klasyczna metoda otwarta Rungego-Kutty czwartego rzędu:
w1 = w4 = 1/6, w2 = w3 = 1/3
k1,n = f(tn, yn), k2,n = f(tn + hn/2, y + hnk1,n/2)
k3,n = f(tn + hn/2, y + hnk2/2), k4,n = f(tn + hn, yn + hnk3)
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 14
Przykład 4
Dla układu RC z Przykładu 1 z liniowym konden-
u(t) dla h=0.25
1
satorem równanie dynamiki ma postać:
0.8
du e(t) - u(t)
metoda Eulera otwarta
= f(t, u), gdzie f(t, u) =
0.6
metoda Eulera zamknieta
dt RC
metoda Heuna
0.4
metoda RK4
Zastosowanie formuły otwartej Eulera prowadzi
rozw. dokladne
0.2
tu do równania różnicowego:
0
0 1 2 3 4 5
t
un+1 = un+h(1-un), gdzie  = 1/(RC)
´[ u(t)]
0
10
Dla formuły zamkniętej Eulera uzyskuje się rów-
nanie różnicowe
-2
10
un+1 = un + h(1 - un+1)
a dla metody Heuna:
-4
10
0 1 2 3 4 5
t
un+1 = un + h(1 - h/2)(1 - un)
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 15
11.3.2 Metody wielokrokowe
Szeroką klasę metod wielokrokowych otrzymuje się przez przedstawienie rozwiązania równania
różniczkowego za pomocą całki:
tn+1
y(tn+1) = y(tn-M) + f(t, y(t))dt
tn-M
a następnie zastąpienie funkcji podcałkowej wielomianem interpolacyjnym opartym na węzłach
wybranych ze zbioru: {tk, f(tk, y(tk))} dla k = n - L, . . . , n + 1.
Uwaga: maksymalny rząd K-krokowej metody zbieżnej wynosi K, gdy K jest liczbą parzystą,
albo K + 1, gdy K jest liczbÄ… nieparzystÄ….
Przykład 5
fn+1-fn
Niech L = 1, M = 0, a przybliżenie liniowe: f(t, y(t)) H" fn + (t - tn)
h
tn+1
fn+1 - fn fn+1 + fn
yn+1 = yn + [fn + (t - tn)]dt = yn + h
h 2
tn
Jest to tzw. formuła trapezów. Formuła jest zamknięta (interpolacyjna), rzędu 2.
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 16
Wykonanie każdego kroku formuły zamkniętej wymaga rozwiązania układu równań
nieliniowych. Przykładowo, dla formuły trapezów jest to układ równań:
f(tn+1, yn+1) + f(tn, yn)
yn+1 = yn + h
2
Dla iteracyjnego rozwiązania tego układu względem yn+1 (np. za pomocą metody
Newtona-Raphsona) potrzeba dobrego przybliżenia początkowego. Można go uzyskać za
pomocą otwartej formuły całkowania tego samego rzędu co dana formuła zamknięta.
Przykład 6
fn-fn-1
Niech L = 1, M = 1, a przybliżenie liniowe: f(t, y(t)) H" fn-1 + (t - tn-1)
h
tn+1
fn - fn-1 3fn - fn-1
yn+1 = yn + [fn-1 + (t - tn-1)]dt = yn + h
h 2
tn
Formuła jest otwarta (ekstrapolacyjna), rzędu 2. Może służyć jako predyktor dla zamkniętej
metody rzędu 2 (np. trapezów), nazywanej wówczas korektorem.
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 17
11.4 Problem do samodzielnej analizy
Dla nieliniowego układu RC pokazanego na rysunku przyjąć, R=1M&!, e(t) = 0 dla t < 0 i
e(t) = 1 dla t 0, q(t) = 0 dla t 0, a także że równanie ładunkowe kondensatora ma
postać:
R
u
Q(u) = C(x)dx, gdzie
0
e(t) q(t) u(t)
C0
C(u) = , m = 0.5, U0 = 0.6, C0 = 10-12
u
(1 + )m
U0
" Zapisać równania różnicowe metody otwartej i metody zamkniętej Eulera dla tego układu
elektrycznego, przyjmujÄ…c zerowe warunki poczÄ…tkowe.
" Rozwiązać w środowiskuMATLABrównania różnicowe otwartej i zamkniętej metody
Eulera. Dla obydwóch metod wykreślić przebieg u(t); porównać z rozwiązaniem
uzyskanym za pomocÄ… metodyode45.
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 18
Wybrane formuły predykcyjno-korekcyjne
Metody Adamsa-Bashfortha p Ap+1
yn+1 = yn + hfn (m. otwarta Eulera) 1 -1/2
yn+1 = yn + h(3fn - fn-1)/2 2 -5/12
yn+1 = yn + h(23fn - 16fn-1 + 5fn-2)/12 3 -3/8
yn+1 = yn + h(55fn - 59fn-1 + 37fn-2 - 9fn-3)/24 4 -251/720
Metody Adamsa-Moultona p Ap+1
yn+1 = yn + hfn+1 (m. zamknięta Eulera) 1 1/2
yn+1 = yn + h(fn+1 + fn)/2 (m. trapezów) 2 1/12
yn+1 = yn + h(5fn+1 + 8fn - fn-1)/12 3 1/24
yn+1 = yn + h(9fn+1 + 19fn - 5fn-1 + fn-2)/24 4 19/720
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 19
11.5 Zbieżność i dokładność dyskretnych metod całkowania
Użyteczna N-krokowa stałokrokowa formuła całkowania równań różniczkowych musi być
zbieżna, tzn. dla każdego warunku początkowego y0 = y(t0)
" błąd globalny a" maxi=0,...N(h) ei dąży do 0 dla kroku całkowania h 0, czyli
g
" limh0 yn = y(tn)
Ważne jest też zachowanie metody dla małych, ale niezerowych wartości kroku h.
Niech dla K-krokowej metody całkowania
K K
yn+1 = Ä…kyn+1-k + h ²k fn+1-k, n = 0, 1, . . . , N - 1
k=1 k=0
y n+1-k
znane są dokładne wartości yi = y(ti), fi = y (ti), i = n, n - 1, n - 1 + K.
Lokalny błąd całkowania w n + 1-tym kroku h można zdefiniować jako:
rn+1 = yn+1 - y(tn+1)
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 20
Jeśli zależność błędu lokalnego od kroku h jest dostatecznie gładka, to:
h2 h3
rn+1 a" rn+1(h) = rn+1(0) + r n+1(0)h + r n+1(0) + r n+1(0) + · · ·
2 3!
Zależność rn+1(h) jest wykorzystywana do określania rzędu metody.
Metoda jest rzędu p, jeśli pierwsze p pochodnych funkcji rn+1(h) jest równe 0. Błąd lokalny
można wówczas wyrazić jako:
r(p+1)(0)
rn+1(h) = y(p+1)hp+1 + O(hp+2)
(p + 1)!
Ap+1
Przykład 7. Wyznaczmy błąd lokalny i rząd otwartej metody Eulera:yn+1 = yn + hf(tn, yn).
Jeśli yn = y(tn), to rn+1 = yn+1 - y(tn+1) = y(tn) + hy (tn) - y(tn + h).
RozwijajÄ…c rn(h) w szereg Taylora:
h2 h2
rn+1(h) = y(tn) + hy (tn) - (y(tn) + y (tn)h + y (tn) + · · ·) = -y (tn) + · · ·
2 2
2
Stąd rząd metody p = 1, a błąd lokalny rn+1(h) H" -y (tn)h
2
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 21
Przykład 8
Przy rozwiązywaniu skalarnego równania różniczkowego za pomocą otwartej metody Eulera,
niedokładność przybliżenia w pierwszym kroku y1 jest obarczone lokalnym błędem obcięcia):
2
y1 - y(t1) = r1 H" -y (tn)h
2
Niedokładność przybliżenia uzyskanego w drugim kroku e2 = y2 - y(t2) nie jest prostą sumą
lokalnego błędu obcięcia r2 oraz błędu globalnego (ozn. e1 = y1 - y(t1)) popełnionego przy
wyznaczaniu y1:
e2 = y2 - y(t2) = y1 + hf(t1, y1) - y(t1 + h) =
h2
= y(t1) + e1 + hf(t1, y(t1) + e1) - (y(t1) + y (t1)h - y (t1) + · · ·) =
2
h2
= e1 + h(f(t1, y(t1) + e1) - y (t1)) + y (t1) + O(h3)
2
h2 h2
|e2| |e1| + hL|e1| + |y (t1)| + O(h3) (1 + hL)|e1| + M2 + O(h3)
2 2
gdzie L jest stałą Lipschitza, a M2 = sup[t0,tF ] y (t).
Za pomocą indukcji zupełnej można udowodnić, że po N krokach
hM2 hM2 L
|eN | ((ehL)N - 1) = (eT - 1) = const · h1), T = Nh
2L 2L
Metoda ta jest wiÄ™c zbieżna (|eN | “! 0 dla h “! 0), rzÄ™du p=1 (gdyż |eN | = O(hp), p = 1).
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 22
Przykład 9
Dobrać współczynniki Ä…, ² tak, by jednokrokowa otwarta formuÅ‚a
yn+1 = Ä…yn + h²f(tn, yn)
była możliwie wysokiego rzędu.
rn+1(h) = yn+1 - y(tn + h) = Ä…y(tn) + h²y (tn) - y(tn + h) =
h2
= Ä…y(tn) + h²y (tn) - (y(tn) + hy (tn) + y (tn) + O(h3)) =
2
h2
= y(tn)(Ä… - 1) + y (tn)h(² - 1) - y (tn) + O(h3)
2
Dla Ä… = 1 i ² = 1 jedynie pierwsza pochodna funkcji rn+1(h) jest równa 0, a wiÄ™c uzyskana
formuła (Eulera) jest rzędu p = 1. Błąd lokalny wyraża się następująco:
h2
rn+1(h) = -y (tn) + O(h3)
2
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 23
Automatyczny wybór kroku całkowania
Załóżmy, że metoda całkowania jest rzędu p i że dominuje lokalny błąd dyskretyzacji rzędu
hp+1. Zgodnie z zasadą Rungego, następny punkt rozwiązania wyznacza się dwukrotnie:
y
n+1,1
"n+1
" yn+1,1 : przy użyciu formuły z krokiem h
y
n+1,2
" yn+1,2: przez dwukrotne użycie formuły z kro-
kiem h/2
t
y
n
t +h
t t +h/2
n n n
yn+1,1 = y(t + h) + cp+1hp+1 + O(hp+2), cp+1 = Ap+1y(p+1)
yn+1,2 = y(t + h) + 2cp+1(h/2)p+1 + O(hp+2)
"n+1 = yn+1,1 - yn+1,2 może służyć do oszacowania lokalnego błędu obcięcia.
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 24
"n+1 = cp+1hp+1 - 2cp+1(h/2)p+1 + O(hp+2) H" cp+1hp+1(1 - 2-p)
"n+1
cp+1 H"
hp+1(1-2-p)
Oszacowanie błędu popełnianego przy obliczaniu yn+1,1 i yn+1,2:
¾n+1,1 = yn+1,1-y(tn+1) H" "n+1/(2p-1), ¾n+1,2 = yn+1,2-y(tn+1) H" 2-p"n+1/(2p-1)
Przykładowy algorytm doboru długości kroku (z podwajaniem)
1. Wyznacz ¾n+1,1, ¾n+1,2 dla dÅ‚ugoÅ›ci kroku hn+1 i hn+1/2
2. JeÅ›li |¾n+1,1| < emax, to krok h należy podwoić
3. JeÅ›li |¾n+1,1| > emax, ale ¾n+1,2 < emax, to należy kontynuować z niezmnienionym
krokiem.
4. JeÅ›li |¾n+1,2| > emax,to należy powtórzyć od p. 1 ze skróconym dwukrotnie krokiem.
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 25
Przykład 10
Wyznaczyć odpowiedz
układu z Przykładu 1 dla R = 1, C = 1, E0 = 1, t0 = 0, T = 5. Zastosować
otwartą metodę Eulera: a) z ustalonym krokiem, b) ze zmiennym krokiem. Założyć, że błąd globalny
|eN| emax = 0.005
Dokładne rozwiązanie: y(t) = 1 - exp(-t), dla t 0
Porownanie bledow calkowania
0
L = supt"[0,T ]{|y(1)(t)|} = 1
10
M2 = maxt"[0,T ]{|y(2)(t)|} = 1
Błąd globalny metody stałokrokowej: -5
10
zmienny krok
exp(5) - 1
staly krok
|eN| h H" 73, 7h
-10
2
10
0 1 2 3 4 5
t [s]
StÄ…d h H" 0, 005/73, 7 = 6, 78 · 10-5
Uzywane dugosci kroku
0
10
Liczba kroków czasowych:
" 73707: algorytm ze stałym krokiem, maks. zmienny krok
staly krok
bÅ‚Ä…d: 1.2 · 10-5 emax
" 373: algorytm z podwajaniem kroku, maks. -5
10
0 1 2 3 4 5
t [s]
bÅ‚Ä…d: 1.37 · 10-3 < emax
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
h [s]
Rozwiązywanie równań różniczkowych 26
11.6 Stabilność dyskretnych metod całkowania
Przykład 11
Wyznaczyć przebiegi czasowe napięć u1(t) i u2(t) dla t " [0; 0.015] i e(t) = 1(t)
(włączenie z
ródła jednostkowego w chwili 0) korzystając z formuły otwartej Eulera. Porównać
wyniki uzyskane dla N = 7, 10, 30 kroków czasowych.
u1(t) u2(t)
4 0.015
N=7
N=10
N=30
przebieg dokladny
R R
3
1 2
2 0.01
e(t)
u1 u2
C C
1
1 2
0 0.005
R1 = 1k&!, R2 = 1M&!,
-1
C1 = C2 = 1µ F,
-2 0
0 0.005 0.01 0.015 0 0.005 0.01 0.015
t t
Dla zbyt długich kroków czasowych w odpowiedzi pojawiają się pasożytnicze oscylacje . Jest to
objaw niestabilności numerycznej algorytmu całkowania. Formuła zamknięta Eulera nie
wykazuje takiego zachowania.
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 27
Analizę stabilności metod całkowania równań różniczkowych wystarczy przeprowadzić dla zadania
testowego: y (t) =  y(t), y(t0) = y0 , gdzie  jest (ogólnie: zespolonym) współczynnikiem.
Dla metod wielokrokowych
K K
yn+1 = Ä…kyn+1-k + h ²kfn+1-k, n = 0, 1, . . . , N - 1
k=1 k=0
użycie równania testowego oznacza przyjęcie: fi = yi, a stąd:
K K K
yn+1 = Ä…kyn+1-k + h ²kyn+1-k = h²0yn+1 + (Ä…k + hyn+1-k)
k=1 k=0 k=1
K
1
yn+1 = (Ä…k + h²k)yn+1-k
1 - h²0 k=1
Rozwiązanie ogólne tego równania różnicowego zawiera składową postaci yn = zn, gdzie z jest liczbą
zespoloną, która jest pierwiastkiem następującego wielomianu charakterystycznego
K
1
Á(z) = zK - (Ä…k + h²k)zK-k
1 - h²0 k=1
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 28
Metoda dyskretna rozwiązywania równań różniczkowych jest stabilna, jeśli generowane przez
niÄ… rozwiÄ…zanie zadania testowego yn sÄ… ograniczone.
Metodę różnicową nazywany absolutnie stabilną w pewnym obszarze &! płaszczyzny zespolonej
C, jeśli wszystkie pierwiastki wielomianu charakterystycznego spełniają warunek |z| 1 (a
pierwiastki wielokrotne: |z| < 1) dla każdego h " &!. Jeśli &! zawiera lewą półpłaszczyznę
płaszczyzny C, to metodę nazywamy A-stabilną.
Przykład 12
Wykonajmy analizę stabilności otwartej metody Eulera: yn+1 = yn + hf(tn, yn).
Dla równania testowego f(tn, yn) =  yn uzyskujemy równanie różnicowe:
yn+1 = yn + hyn, czyli yn+1 = yn(1 + h)
Rozwiązanie ogólne tego równania różnicowego zawiera składową postaci yn = zn, gdzie z
jest liczbą zespoloną. Stąd z spełnia następujące równanie charakterystyczne: z = 1 + h.
Dla otwartej metody Eulera mamy więc warunek stabilności: |1 + h| 1
Uwaga. Wydłużając odpowiednio krok h można zawsze spowodować niespełnienie tego
warunku nawet dla stabilnego równania testowego (tzn. gdy Re() < 0).
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 29
Przykład 13
Wykonajmy analizę stabilności zamkniętej metody Eulera:
yn+1 = yn + hf(tn+1, yn+1)
Dla równania testowego f(tn, yn) =  yn uzyskujemy równanie różnicowe:
yn+1 = yn + hyn+1, czyli yn+1 = yn/(1 - h)
Równanie charakterystyczne ma tu postać: z = 1/(1 - h), a warunek stabilności:
1/|1 - h| 1
Dla stabilnego równania testowego Re() < 0, co oznacza, że warunek stabilności jest
spełniony dla każdego h 0. W tej sytuacji krok całkowania można wybierać możliwie długi
(ze względu na nakłady obliczeniowe), ograniczony jedynie ze względu na pożądaną
dokładność rozwiązania. Dla metody otwartej Eulera krok jest ograniczony również przez
warunki stabilności.
Podobne relacje, jak dla otwartej i zamkniętej metody Eulera, zachodzą dla metod otwartych i
zamkniętych innych rzędów.
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 30
Przykład 14
Brzegi obszarów stabilności K-krokowych formuł predykcyjnych Adamsa-Bashfortha i formuł
korekcyjnych Adamsa-Moultona pokazano na wykresach.
Algorytmy Adamsa-Bashfortha Algorytmy Adamsa-Moultona
1 5
0.8 4
0.6 3 K=2
K=3
0.4 2
K=4
K=1
0.2 1
Im ( h ) K=1 K=2 K=3 K=4 Im ( h )
0 0
-0.2 -1
-0.4 -2
-0.6 -3
-0.8 -4
-1 -5
-2 -1.5 -1 -0.5 0 0.5 -8 -6 -4 -2 0 2
Re ( h ) Re ( h )
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 31
11.7 Równania sztywne
Układ liniowych równań różniczkowych y = Ay nazywany jest sztywnym, jeśli
n
1
1
gdzie 1 < 2 < · · · < n sÄ… wartoÅ›ciami wÅ‚asnymi macierzy A.
Przykład 15. Równania układu R2C2 (Przykład 3)
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
3
ðÅ‚u1(t)ûÅ‚=ðÅ‚-1001 1 ûÅ‚ ðÅ‚u1(t)ûÅ‚ ðÅ‚10 · 1(t)ûÅ‚
+
u 2(t) 1 -1 u2(t) 0
A
są sztywne, gdyż
"wartości własne macierzy A: "
1 = -501 + 250001 H" -0.999, 2 = -501 - 250001 H" -1001 sÄ… znacznie
różniące się od siebie.
Interpretacja elektryczna sztywności układu: w odpowiedziach czasowych istnieją składowe
bardzo wolne i szybkie.
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 32
Całkowanie sztywnych układów równań różniczkowych jest trudne za pomocą metod, które mają
const
górne ograniczenia stabilnościowe na wielkość kroku całkowania, np. o postaci: h < dla
maxj|j|
j " R (odpowiedzi aperiodyczne).
Maksymalna długość kroku, wymagana przez stabilność, może być dużo mniejsza od długości kroku
wynikającej z wymagań na dokładność rozwiązania.
Przedział całkowania jest zazwyczaj wyznaczony przez ustalenie się najwolniejszej składowej
rozwiÄ…zania, np. T = 5 maxi Äi = 5/ mini |i|
Obliczanie odpowiedzi sztywnych układów równań różniczkowych może wymagać dużej liczby
kroków caÅ‚kowania Ò! bardzo dużych nakÅ‚adów obliczeniowych.
Przykład 16. Dla układu równań z Przykładu 3 minimalna długość kroku metody otwartej Eulera wynika
z warunku stabilności: |1 + hi| 1, przy czym 1 H" -0.999, 2 H" -1001. Stąd maksymalna
dÅ‚ugość kroku caÅ‚kowania dla metody staÅ‚okrokowej: hmax = mini -2 H" 1.998 · 10-3 jest
i
wyznaczona przez krótszÄ… z dwóch staÅ‚ych czasowych ukÅ‚adu, t.j. Ä1 = 1/|1| H" 1.001,
Ä2 = 1/|2| H" -0.999 · 10-3. Czas caÅ‚kowania, konieczny do uzyskania stanu ustalonego w
ukÅ‚adzie, jest z kolei wyznaczony przez dÅ‚uższÄ… ze staÅ‚ych czasowych: T H" 5 maxi Äi = 5Ä1 H" 5.
Wyznaczenie odpowiedzi wymaga wiÄ™c co najmniej N = T/hmax = 5Ä1/Ä2 = 52/1 H" 2503
kroków całkowania.
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 33
Metody A-stabilne, to metody absolutnie stabilne dla obszaru stabilności w obszarze:
&! = {z " C | Re(z) < 0}
A-stabilność oznacza więc brak ograniczeń kroku całkowania dla rozwiązywania stabilnych równań
liniowych (wartości własne w lewej półpłaszczyz
nie, tzn. dla  " &!.
25
Niestety maksymalny rzÄ…d A-stabilnych (nie posiadajÄ…cych sta-
20
bilnościowych ograniczeń długości kroku) metod wielokrokowych
15
10
Adamsa-Moultona wynosi 2.
5
Istnieją jednak inne metody rzędu >2, np. zamknięte metody Geara,
Im(h)
Obszar
1 2 3 K=4 K=5 K=6
0
stabil-
nosci
-5
których obszar stabilności jest nieograniczonym z lewej strony
-10
podobszarem lewej półpłaszczyzny (w układzie współrzędnych
-15
Re(h), Im(h)).
-20
-25
-10 0 10 20 30
Re(h)
Metody Geara bywają wykorzystywane w symulatorach układów elektronicznych (np.SPICE) do
rozwiązywania sztywnych układów równań o odpowiedziach aperiodycznych oraz umiarkowanie
oscylacyjnych (tłumionych).
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 34
Otwarte metody Geara p Ap+1
yn+1 = hfn + yn (m. otwarta Eulera) 1 -1
2
yn+1 = 2hfn + yn-1 2 -1
3
3 1
yn+1 = 3hfn - yn + 3yn-1 - yn-2 3 -1
2 2 4
10 1
yn+1 = 4hfn - yn + 6yn-1 - 2yn-2 + yn-3 4 -1
3 3 5
65 5 1
yn+1 = 5hfn - yn + 10yn-1 - 5yn-2 + yn-3 - yn-4 5 -1
12 3 4 6
77 3 1
yn+1 = 6hfn - yn + 15yn-1 - 10yn-2 + 5yn-3 - yn-4 + yn-5 6 -1
10 2 5 7
Zamknięte metody Geara p Ap+1
1
yn+1 = hfn+1yn (m. zamknięta Eulera) 1
2
2 4 1 2
yn+1 = hfn+1 + yn - yn-1 2
3 3 3 9
6 18 9 2 3
yn+1 = hfn+1 + yn - yn-1 + yn-2 3
11 11 11 11 22
12 48 36 16 3 12
yn+1 = hfn+1 + yn - yn-1 + yn-2 - yn-3 4
25 25 25 25 25 125
60 300 300 200 75 12 10
yn+1 = hfn+1 + yn - yn-1 + yn-2 - yn-3 + yn-4 5
137 137 137 137 137 137 137
60 360 450 400 225 72 10 20
yn+1 = hfn+1 + yn - yn-1 + yn-2 - yn-3 + yn-4 - yn-5 6
147 147 147 147 147 147 147 343
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008
Rozwiązywanie równań różniczkowych 35
11.8 Zakres stosowania algorytmów całkowania programu MATLAB 7
Algorytm Równania sztywne? Dokładność
ode45 nie średnia
ode23 nie niska
ode113 nie niska do wysokiej
ode15s tak niska do średniej
ode23s tak niska
ode23t umiarkowanie sztywne niska
ode23tb tak niska
L.J. Opalski,  Wstęp do metod numerycznych - materiały do wykładu wersja: 29.XII.2008