Metody numeryczne i statystyka
Rozwiązywanie równań
ró\
niczkowych
Wersja robocza
Ró\
niczkowanie numeryczne  metoda
ró\
niczkowania jednopunktowego w przód
Pierwsza pochodna funkcji f (x) zdefiniowana jest jako
Chcemy obliczyć pochodną funkcji f (x).
f (x + h) - f (x)
f '(x) = lim
h0
h
Zagadnienie to mo\
na rozwiązać numerycznie poprzez
proste przybli\
enie przy ustalonej wielkości kroku h
y(x0 + h) - y(x0 )
y'(x0 ) H"
h
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 2
Ró\
niczkowanie numeryczne  metoda
ró\
niczkowania jednopunktowego w przód
Pierwsza pochodna funkcji f (x) zdefiniowana jest jako
Chcemy obliczyć pochodną funkcji f (x).
f (x + h) - f (x)
f '(x) = lim
h0
h
Zagadnienie to mo\
na rozwiązać numerycznie poprzez
proste przybli\
enie przy ustalonej wielkości kroku h
y(x0 + h) - y(x0 )
y'(x0 ) H"
h
W celu analizy popełnianego błędu wyra\
enie y (x0+ h) mo\
na rozwinąć w szereg Taylora
w zale\
ności od x0
h2 h3
y(x0 + h)= y(x0)+ h Å" y'(x0)+ Å" y''(x0 ) + Å" y'''(x0 ) + ...
2! 3!
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 3
Ró\
niczkowanie numeryczne  metoda
ró\
niczkowania jednopunktowego w przód
Pierwsza pochodna funkcji f (x) zdefiniowana jest jako
Chcemy obliczyć pochodną funkcji f (x).
f (x + h) - f (x)
f '(x) = lim
h0
h
Zagadnienie to mo\
na rozwiązać numerycznie poprzez
proste przybli\
enie przy ustalonej wielkości kroku h
y(x0 + h) - y(x0 )
y'(x0 ) H"
h
W celu analizy popełnianego błędu wyra\
enie y (x0+ h) mo\
na rozwinąć w szereg Taylora
w zale\
ności od x0
h2 h3
y(x0 + h)= y(x0)+ h Å" y'(x0)+ Å" y''(x0 ) + Å" y'''(x0) + ...
2! 3!
Po odjęciu od obu stron y (x0) i podzieleniu obu stron przez h otrzymujemy
y(x0 + h)- y(x0) h h2
= y'(x0)+ Å" y''(x0) + Å" y'''(x0) + ...
h 2 3!
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 4
Ró\
niczkowanie numeryczne  metoda
ró\
niczkowania jednopunktowego w przód
Pierwsza pochodna funkcji f (x) zdefiniowana jest jako
Chcemy obliczyć pochodną funkcji f (x).
f (x + h) - f (x)
f '(x) = lim
h0
h
Zagadnienie to mo\
na rozwiązać numerycznie poprzez
proste przybli\
enie przy ustalonej wielkości kroku h
y(x0 + h) - y(x0 )
y'(x0 ) H"
h
W celu analizy popełnianego błędu wyra\
enie y (x0+ h) mo\
na rozwinąć w szereg Taylora
w zale\
ności od x0
h2 h3
y(x0 + h)= y(x0)+ h Å" y'(x0)+ Å" y''(x0 ) + Å" y'''(x0) + ...
2! 3!
Po odjęciu od obu stron y (x0) i podzieleniu obu stron przez h otrzymujemy
y(x0 + h)- y(x0) h h2
= y'(x0)+ Å" y''(x0 ) + Å" y'''(x0 ) + ... = y'(x0)+ O(h) , gdzie O(h)  bÅ‚Ä…d obciÄ™cia
h 2 3!
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 5
Ró\
niczkowanie numeryczne  metoda
ró\
niczkowania trójpunktowego
Przybli\
one wartości pochodnej funkcji mo\
na tak\
e obliczyć posługując się algorytmem
ró\
niczkowania trójpunktowego. Mo\
na z niego skorzystać gdy wartości zmiennej
niezale\
nej podawane są ze stałym krokiem h i znane są wartości zmiennej zale\
nej
y(x0 - h), y(x0), y(x0 + h)
Wartość pochodnych określa się w przybli\
eniu jako
1
y'(x0 - h) = [- 3y(x0 - h) + 4y(x0 ) - y(x0 + h)]
2h
1
y'(x0 ) = [- y(x0 - h) + y(x0 + h)]
2h
1
y'(x0 + h) = [y(x0 - h) - 4y(x0 ) + 3y(x0 + h)]
2h
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 6
Ró\
niczkowanie numeryczne  metoda
ró\
niczkowania trójpunktowego
W praktyce korzysta się ze wzorów zapisanych iteracyjnie:
1
'
y1 = [- 3y1 + 4y2 - y3]
- pierwszy punkt
2h
1
yi' = [- yi-1 + yi+1] , i = 2,3,...,n -1
- pozostałe punkty
2h
1
'
yn = [yn-2 - 4yn-1 + 3yn]
- ostatni punkt
2h
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 7
Rozwiązywanie równań ró\
niczkowych
 metoda Eulera
Poszukiwane jest rozwiÄ…zanie
równania ró\
niczkowego postaci:
dy
= f (x, y(x)), x "[x0 ,b]
dx
z warunkiem poczÄ…tkowym
y(x0 ) = y0
RozwiÄ…zaniem omawianego problemu
jest poszukiwana wartość funkcji
y = y(x)
w
x = x1 = x0 + h
Przyjmijmy, \
e y* = y (x1) i krzywa (A0A1*) sÄ… rozwiÄ…zaniem zadania. ZnajÄ…c punkt A0(x0, y0)
(warunek poczÄ…tkowy) mo\
na obliczyć wartość funkcji f w tym punkcie. Jest to współczynnik
kierunkowy stycznej do y (x) w punkcie x0. Następnie wyznaczany jest punkt przecięcia tej
stycznej z prostą x = x1 i oznaczany przez A1(x1, y1). Takie postępowanie prowadzi do
powstania błędu metody przy wyznaczaniu wartości yi
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 8
Rozwiązywanie równań ró\
niczkowych
 metoda Eulera
Spróbujmy sformalizować sposób obliczania yi
w postaci wzoru. Z trójkąta (A0A1P) mamy:
y1 - y0
= f (x0 , y0 )
h
stÄ…d
y1 = y0 + h Å" f (x0 , y0 )
Chcąc wyznaczyć wartość funkcji y2
w punkcie x2 = x1 + h = x0 + 2 h
postępujemy analogicznie, tylko zamiast
korzystać z warunku początkowego korzystamy
z uprzednio obliczonych wartości (x1, y1).
Uogólniając takie postępowanie mo\
na zapisać:
yi = yi-1 + h Å" f (xi-1, yi-1), i = 1,2,...
gdzie yi jest wartością funkcji stanowiącej rozwiązanie równania ró\
niczkowego w xi
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 9
Metoda Eulera
MetodÄ™ Eulera mo\
na otrzymać w inny sposób. Zastępując funkcję y = y (x)
w x = x0 + h jej rozwinięciem w szereg Taylora otrzymujemy:
h2
y(x0 + h)= y(x0)+ h Å" y'(x0)+ Å" y''(x0 ) + ...
2!
Wykorzystując dwa pierwsze wyrazy rozwinięcia w szereg Taylora:
y(x0 + h) E" y(x0 )+ h Å" y'(x0 )
co mo\
na inaczej napisać:
y1 = y0 + h Å" f (x0 , y0 )
Podany wzór jest identyczny jak uzyskany uprzednio. Powtarzając takie postępowanie jak
w poprzednim wyprowadzeniu wzór ten mo\
na uogólnić
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 10
Metoda Rungego-Kutty rzędu II
Przypuśćmy, \
e znamy y (xi -1) i chcemy wyznaczyć przybli\
enie yi wartości y (xi -1+ h) . Idea
metod Rungego-Kutty polega na obliczeniu wartości f (x, y) w pewnych szczególnie
dobranych punktach le\
Ä…cych w pobli\
u krzywej rozwiÄ…zania w przedziale (xi -1, xi -1+ h) oraz
utworzeniu takiej kombinacji tych wartości, która z dobrą dokładnością daje przyrost yi  yi -1
Rozwa\
my metodÄ™ Rungego-Kutty II
rzędu, opartą na przybli\
eniu rozwiÄ…zania
równania ró\
niczkowego za pomocÄ…
krzywej opisanej wielomianem stopnia II
(parabolÄ…).
Równanie ró\
niczkowe ma postać
dy
= f (x, y(x)) , y(x0) = y0
dx
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 12
Metoda Rungego-Kutty rzędu II
Korzystając z własności paraboli, dla której współczynnik kierunkowy
siecznej T równy średniej arytmetycznej współczynników
kierunkowych stycznych i w punktach o współrzędnych odpowiednio
(x0, y0), (x1, y1).
Pozwala to na przybli\
enie rozwiązania równania ró\
niczkowego
łamaną utworzoną z siecznych, a nie stycznych jak miało to miejsce
w metodzie Eulera.
Poniewa\
nie jest znana postać rozwiązania rzędną y1 przybli\
a siÄ™
E
y1
wartością wyznaczoną z metody Eulera
Mo\
na zapisać przybli\
enie poszukiwanego rozwiÄ…zania w punkcie
o współrzędnych (x1, y1) w sposób następujący
S0 + S1
öÅ‚
y1 = y0 + hëÅ‚
ìÅ‚ ÷Å‚
2
íÅ‚ Å‚Å‚
Podstawiając za współczynnik kierunkowy stycznej S0 wyra\
enie
f (x0, y0 ) oraz za współczynnik kierunkowy stycznej S1 wyra\
enie
E
y1
f (x1, ) otrzymuje siÄ™
1
E
y1 = y0 + h[f (x0, y0)+ f (x1, y1 )]
2
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 13
Metoda Rungego-Kutty rzędu II
Pamiętając o tym, \
e:
E
x1 = x0 + h oraz y1 = hf (x0, y0)
wyra\
enie na y1 mo\
na zapisać następująco
1
y1 = y0 + [h Å" f (x0, y0)+ hÅ" f (x0 + h, y0 + hÅ" f (x0, y0))]
2
OznaczajÄ…c
k1 = hf (x0, y0)
Wyra\
enie przyjmuje postać
1
y1 = y0 + [k1 + hÅ" f (x0 + h, y0 + k1)]
2
h Å" f (x0 + h, y0 + k1)
Zastępując wyra\
enie współczynnikiem k2 otrzymuje się wzór na
metodÄ™ Rungego-Kutty II w postaci
1
y1 = y0 + [k1 + k2]
2
k1 = hf (x0, y0)
k2 = hf (x0 + h, y0 + k1)
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 14
Metoda Rungego-Kutty rzędu II
Wyznaczając kolejne wartości poszukujemy funkcji yi stanowiącej przybli\
one rozwiÄ…zanie
równania ró\
niczkowego w oparciu o sieczne według przedstawionej metodyki wzór metody
Rungego-Kutty II przyjmuje ogólną postać:
1
yi = yi-1 + [k1 + k2]
2
k1 = hf (xi-1, yi-1)
k2 = hf (xi-1 + h, yi-1 + k1)
dla i = 1,2,..., N
Ogólnie metody Rungego-Kutty polegajÄ… na takim dobraniu współczynników a, b, & , Ä…, ², Å‚,&
oraz liczb R1, R2,& tak aby wartość yi określona przez ciąg równań była mo\
liwie jak najbli\
sza
dokładnej wartości:
yi = yi-1 + R1k1 + R2k2 +...+ Rnkn , i =1,2..., n =1,2,...
k1 = hf (xi-1, yi-1)
k2 = hf (xi-1 + ah, yi-1 +Ä…k1)
k3 = hf (xi-1 + bh, yi-1 + ²k1 + Å‚k2)
...
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 15
Metoda Rungego-Kutty rzędu IV
Najbardziej znana jest metoda Rungego-Kutty czwartego rzędu opisana zale\
nością
1
yi = yi-1 + (k1 + 2 Å" k2 + 2 Å" k3 + k4 )
6
gdzie
k1 = h Å" f (xi-1, yi-1)
k1
h
k2 = h Å" f (xi-1 + , yi-1 + )
2 2
h k2
k3 = h Å" f (xi-1 + , yi-1 + )
2 2
k4 = h Å" f (xi-1 + h, yi-1 + k3 )
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 16
Rozwiązywanie równań ró\
niczkowych
pierwszego rzędu
Przykład 1 - Rozładowanie kondensatora przez opornik
Rozpatrzmy obwód elektryczny w którym kondensator
W
o pojemności C został naładowany do napięcia U =U0.
W chwili t = 0 został zwarty włącznik W i nastąpiło
rozładowanie kondensatora przez oporność R
Q dQ dU
U = , I = = C
, Q - Å‚adunek
C dt dt
dU U
C + = 0
dt R
Jest to równanie ró\
niczkowe zwyczajne rzędu I . Aby je rozwiązać dzielimy obie
strony równania przez C oraz po lewej stronie pozostawiamy najwy\
szÄ… pochodnÄ… U
dU U
= -
dt RC
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 17
Algorytm metody Eulera
Etap:
Start
1. Zdefiniowanie funkcji f (x, y (x)), przyjęcie granic
f (x, y (x)),
x0, xk, N, y0 badanego przedziału [x0, xk], liczby kroków N, wartości
poczÄ…tkowej y0
h = (xk - x0) / N
2. Obliczenie długości kroku h, dyskretyzacja przedziału
xi=h·i, i=1,2,...,N
[x0, xk]
nie
id"N
tak
3. Pętla obliczająca rekurencyjnie i-te wartości yi na
yi+1=yi+h·f(xi,y(xi))
podstawie wzoru Eulera
yi = yi-1 + h Å" f (xi-1, yi-1), i = 1,2,..., N
i=i+1
y(xi)
4. Wynik  dyskretna funkcja y (xi ) będąca przybli\
eniem
rozwiązania równania ró\
niczkowego
Koniec
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 18
Rozwiązywanie równań ró\
niczkowych
pierwszego rzędu
Rozwiązanie metodą Eulera dla danych: napięcie początkowe U0 = 10 [V] w chwili t0 = 0 [s],
oporność R = 10 [&!], pojemność C = 0,1 [F], czas końcowy tk = 2 [s], ilość kroków N = 4
Dla równania ró\
niczkowego o postaci
dy
= f (x, y(x)), x "[x0 , xk ]
dx
Wzór iteracyjny Eulera ma postać:
yi = yi-1 + h Å" f (xi-1, yi-1), i = 1,2,..., N
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 19
Rozwiązywanie równań ró\
niczkowych
pierwszego rzędu
Rozwiązanie metodą Eulera dla danych: napięcie początkowe U0 = 10 [V] w chwili t0 = 0 [s],
oporność R = 10 [&!], pojemność C = 0,1 [F], czas końcowy tk = 2 [s], ilość kroków N = 4
Dla równania ró\
niczkowego o postaci
dy
= f (x, y(x)), x "[x0 , xk ]
dx
Wzór iteracyjny Eulera ma postać:
yi = yi-1 + h Å" f (xi-1, yi-1), i = 1,2,..., N
Za funkcjÄ™ f (x, y(x)) mo\
na podstawić:
1
a = - = -1 , y(x) = U (t)
, gdzie
f (x, y(x)) = a Å" y(x)
RC
W przedziale [x0,xk] rozwiązanie jest badane dla wartości xi = xi-1+ h, gdzie x0 = 0, xk = 2,
h = 0,5
Wartość początkowa y0 = U0 = 10
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 20
Rozwiązywanie równań ró\
niczkowych
pierwszego rzędu
Dla pierwszego kroku obliczeń (i=1):
y1 = y0 + h Å" a Å" y0 = 0 - 0,25 Å"1Å"10 = 10
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 21
Rozwiązywanie równań ró\
niczkowych
pierwszego rzędu
Dla pierwszego kroku obliczeń (i=1):
y1 = y0 + h Å" a Å" y0 = 0 - 0,25 Å"1Å"10 = 10
Dla drugiego kroku obliczeń (i=2):
y2 = y1 + h Å" a Å" y1 = 10 - 0,5 Å"1Å"10 = 5
Dla trzeciego kroku obliczeń (i=3):
y3 = y2 + h Å" a Å" y2 = 5 - 0,25Å"1Å" 5 = 0,25
Dla czwartego kroku obliczeń (i=4):
y4 = y3 + h Å" a Å" y3 = 0,25 - 0,25 Å"1Å"10 = 0.625
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 22
Rozwiązywanie równań ró\
niczkowych
pierwszego rzędu
Dla pierwszego kroku obliczeń (i=1):
y1 = y0 + h Å" a Å" y0 = 0 - 0,25 Å"1Å"10 = 10
Dla drugiego kroku obliczeń (i=2):
y2 = y1 + h Å" a Å" y1 = 10 - 0,5 Å"1Å"10 = 5
Dla trzeciego kroku obliczeń (i=3):
y3 = y2 + h Å" a Å" y2 = 5 - 0,25Å"1Å" 5 = 0,25
Dla czwartego kroku obliczeń (i=4):
y4 = y3 + h Å" a Å" y3 = 0,25 - 0,25 Å"1Å"10 = 0.625
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 23
Rozwiązywanie równań ró\
niczkowych
wy\
szych rzędów
Przedstawione wcześniej metody pozwalają na rozwiązywanie równań ró\
niczkowych
zwyczajnych I rzędu. Równanie ró\
niczkowe wy\
szego rzędu mo\
na rozwiązać po
sprowadzaniu go do układu równań ró\
niczkowych rzędu I
W przypadku rozwiązywania równań ró\
niczkowych wy\
szego rzędu
(
y(n) = f (x, y, y',..., y(n-1)), y(x0 ) = y0 , y'(x0 ) = y'0 ,..., y(n-1) (x0 ) = y0n-1)
równanie zostaje sprowadzone do układu równań rzędu I przy wykorzystaniu metody
podstawienia zmiennych. StosujÄ…c notacjÄ™
y1 = y, y2 = y', y3 = y'', .... , yn = y(n-1)
równowa\
ny układ równań przybiera postać
y1'= y2
Å„Å‚
ôÅ‚y '= y3
ôÅ‚
2
òÅ‚
...
ôÅ‚
ôÅ‚yn '= f y1, y2 ,..., yn
(x, )
ół
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 24
Rozwiązywanie równań ró\
niczkowych
wy\
szych rzędów
Przykład 2
Rozpatrzmy układ drgający o jednym stopniu swobody x (t), masie m, współczynniku
sprÄ™\
ystości k, współczynniku tłumienia b. Wymuszeniem jest siła harmoniczna
F(t)=F·sin(Ét) , gdzie F oznacza amplitudÄ™, É  czÄ™stość drgaÅ„
Równanie ruchu opisujące ten układ ma postać:
&& &
mx + bx + kx = F(t)
Jest to równanie ró\
niczkowe zwyczajne rzędu II. Aby je
sprowadzić do układu dwóch równań ró\
niczkowych rzędu I
dzielimy je stronami przez m
b k F(t)
&& &
x + x + x =
m m m
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 25
Rozwiązywanie równań ró\
niczkowych
wy\
szych rzędów
Po lewej stronie pozostawiamy najwy\
szÄ… pochodnÄ… x
b k F(t)
&& - x - x +
&
x =
m m m
Dokonujemy następującego podstawienia:
x = x1
Å„Å‚
òÅ‚x = x2
&1
ół
Równanie ró\
niczkowe mo\
na zapisać jako układ dwóch równań rzędu I:
&
x1 = x2
Å„Å‚
ôÅ‚
òÅ‚ k b F(t)
&2
ôÅ‚x = - m x1 - m x2 + m
ół
Układ równań zapisany w postaci macierzowej:
0 1
îÅ‚ Å‚Å‚
&
x1 x1 îÅ‚ 0 Å‚Å‚
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
ïÅ‚ śł ïÅ‚
= Å" +
F(t)śł
ïÅ‚x śł ïÅ‚x śł
ïÅ‚- k - b śł ïÅ‚ śł
&2
ðÅ‚ ûÅ‚ ðÅ‚ 2 ûÅ‚
ðÅ‚ m m ðÅ‚ m ûÅ‚
ûÅ‚
&
X = A Å" X + B
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 26
Zastosowanie funkcji środowiska MATLAB do
rozwiązywania równań ró\
niczkowych
Dane: Masa m = 1 [kg], sprÄ™\
ystość k = 100, tłumienie b = 0,1. Wymuszenie  siła
harmoniczna F (t) = F·sin(Ét) o amplitudzie F = 1 [N]. Szukane jest przemieszczenie x
masy m w okolicy częstości rezonansowej układu w przedziale czasu [t0 = 0, tk = 10 ] przy
zadanych warunkach poczÄ…tkowych x0 = 0, x 0 = 0
Aby rozwiązać układ przy pomocy funkcji środowiska Matlab nale\
y:
&
X
1. Napisać m-plik funkcyjny zwracający wartości funkcji f (x, t ) zapisanej macierzowo jako ,w którym:
" Przyjęte zostają stałe parametry opisujące układ m, k, F, b
" Na podstawie wyznaczonego wcześniej układu równań ró\
niczkowych (zapisanego w postaci
macierzowej) obliczane są wartości prawej strony równania
0 1 0
îÅ‚ Å‚Å‚ îÅ‚ Å‚Å‚
&
x1
îÅ‚ Å‚Å‚ x1
îÅ‚ Å‚Å‚
,
& ïÅ‚ śł ïÅ‚
X = = A Å" X + B A = , X = , B =
F(t)śł,
ïÅ‚x śł
ïÅ‚x śł
ïÅ‚- k - b śł ïÅ‚ śł
&2
ðÅ‚ 2 ûÅ‚
ðÅ‚ ûÅ‚
ðÅ‚ m m ðÅ‚ ûÅ‚
ûÅ‚ m
2. Napisać m-plik skryptowy definiujący warunki początkowe oraz rozwiązujący równanie
ró\
niczkowe, w którym:
" Określony jest przedział argumentów rozwiązania [t0, tk]
" Wprowadzone sÄ… warunki poczÄ…tkowe x0 , x 0
" Rozwiązane jest układ równań zapisany w pliku funkcyjnym przy wykorzystaniu funkcji
środowiska (ode23, ode45)
" Przedstawione jest rozwiÄ…zanie w postaci wykresu funkcji x(t)
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 27
Zastosowanie funkcji środowiska MATLAB do
rozwiązywania równań ró\
niczkowych
Składnia funkcji (ode23, ode45,& ) rozwiązującej równania ró\
niczkowe zwyczajne (ODE  ang. ordinary
differential equations) rzędu I jest następująca (w MATLAB 6.5):
[X,Y] = solver(odefun,tspan,y0,options)
gdzie solver oznacza nazwÄ™ funkcji reprezentujÄ…cej algorytm rozwiÄ…zania. Funkcje ode23 i ode45
korzystają z par metod Rungego-Kutty odpowiednio: rzędu 2 i 3 (ode23) oraz 4 i 5 (ode45). Funkcje te
zwracają wektor Y stanowiący rozwiązanie równania oraz X  wektor argumentów funkcji, który wyznaczany
jest adaptacyjnie na podstawie zało\
onej tolerancji błędu
odefun  odwołanie pośrednie do funkcji znajdującej się po prawej stronie równania ró\
niczkowego przy
pomocy znaku @, na przykład - @nazwa, gdzie nazwa oznacza nazwę funkcji
tspan - wektor określający granice przedziału w którym rozwiązywane jest równanie [x0 xk]
y0  wektor określający warunki początkowe y(x0)
options  opcjonalny parametr zawierający parametry całkowania, stworzony przy pomocy funkcji odeset
options = odeset('name1',value1,'name2',value2,...)
A
ańcuch tekstowy 'name1', 'name2',& określa parametr, natomiast value1, value2,& jego wartość.
Przykładowo za  name1 ,  name2 ,& mo\
na przyjąć
AbsTol  okreÅ›lajÄ…cy tolerancjÄ™ bÅ‚Ä™du bezwzglÄ™dnego, domyÅ›lna wartość to 1·10-6
RelTol  okreÅ›lajÄ…cy tolerancjÄ™ bÅ‚Ä™du wzglÄ™dnego, domyÅ›lna wartość to 1·10-3
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 28
Zastosowanie funkcji środowiska MATLAB do
rozwiązywania równań ró\
niczkowych
W celu stworzenia pliku funkcyjnego o nazwie NazwaFunkcji.m w oknie poleceń (Command
Window) nale\
y wpisać
>> edit NazwaFunkcji
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 29
Zastosowanie funkcji środowiska MATLAB do
rozwiązywania równań ró\
niczkowych
Otwarte zostaje okno edytora, w którym nale\
y umieścić kod pliku funkcyjnego
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 30
Zastosowanie funkcji środowiska MATLAB do
rozwiązywania równań ró\
niczkowych
m-plik funkcjny  NazwaFunkcji.m opis
Wyra\
enie definiujÄ…ce funkcjÄ™ o nazwie
function Xprim = NazwaFunkcji(t,x)
NazwaFunkcji w której parametry wejściowe to
t, x oraz parametrem wyjściowym jest Xprim
m=1 % masa [kg]
Przyjęcie stałych parametrów równania
F=1 % amplituda siły [N]
k=100 % wsp. sprÄ™\
ystości
b=.2 % wsp. tłumienia
Przyjęcie częstości wymuszenia równej
omega=sqrt(k/m) % częstość [rad/s]
częstości drgań własnych układu
A=[ 0 1
Wprowadzenie macierzy opisujących równanie
-k/m -b/m];
ruchu układu
X=[ x(1)
x(2)];
B=[ 0
F*sin(omega*t)/m];
Xprim=A*X+B; % wynik
Obliczenie Xprim
Kolor zielony  komentarz
Kolor niebieski  wyra\
enie Matlab
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 31
Zastosowanie funkcji środowiska MATLAB do
rozwiązywania równań ró\
niczkowych
Po wprowadzeniu kodu plik nale\
y zapisać w katalogu roboczym File -> Save (lub Ctr+S)
Wa\
ne! Nazwa podana przy definiowaniu funkcji musi być nazwą m-pliku
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 32
Zastosowanie funkcji środowiska MATLAB do
rozwiązywania równań ró\
niczkowych
W celu stworzenia pliku skryptowego o nazwie NazwaSkryptu w oknie poleceń nale\
y wpisać
>> edit NazwaSkryptu
co powoduje, \
e w kolejnej zakładce otwarty jest nowy plik tekstowy o nazwie NazwaSkryptu
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 33
Zastosowanie funkcji środowiska MATLAB do
rozwiązywania równań ró\
niczkowych
m-plik skryptowy  NazwaSkryptu.m opis
t0=0 %czas poczÄ…tkowy
Przyjęcie przedziału czasu oraz parametrów
początkowych równania
tk=10 %czas końcowy
x0=0 %przemieszczenie poczÄ…tkowe
xprim0=0 %prędkość początkowa
Rozwiązanie równania ró\
niczkowego przy
[T,X]=ode45(@NazwaFunkcji,[t0 tk],...
pomocy funkcji ode45, przy ustawieniach
[x0 xprim0]);
domyślnych dla funkcji zapisanej w pliku
NazwaFunkcji.m
figure(1) %nowe okno graficzne
Otworzenie okna graficznego. Narysowanie
plot(T,X(:,1)) %wykres funkcji x(t)
wykresu funkcji przemieszczeń x(t)
xlabel('t - czas [s]')
Opisanie osi wykresu
ylabel('x - przemieszczenie [m]')
Kolor zielony  komentarz
Kolor czerwony  łańcuch tekstowy
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 34
Zastosowanie funkcji środowiska MATLAB do
rozwiązywania równań ró\
niczkowych
Po wprowadzeniu kodu plik nale\
y zapisać w katalogu roboczym File -> Save (lub Ctr+S)
Program uruchamia siÄ™ wybierajÄ…c w menu Debug -> Run (lub klawisz funkcyjny F5)
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 35
Zastosowanie funkcji środowiska MATLAB do
rozwiązywania równań ró\
niczkowych
Efektem działania programu jest wykres przemieszczeń x w zale\
ności od czasu t
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 36
Bezpośrednie całkowanie równań ruchu.
Metoda Newmarka
W roku 1956 Newmark zaproponował rodzinę metod umo\
liwiajÄ…cych jednokrokowe
rozwiązywanie liniowych dynamicznych równań ruchu o postaci
&& &
Mxt + Cxt + Kxt = Ft
&
xt
Przy zało\
eniu obliczeń w chwilach czasu "t, 2"t, 3"t,& rozwinięcie wyra\
enia xt oraz
w szereg Taylora w punkcie t prowadzi do uzyskania równań:
2 3
("t) ("t)
& &&t -"t &&&t -"t
xt = xt -"t + "txt -"t + x + x + ...
2 6
2
("t)
& & && &&&t -"t
xt = xt -"t + "txt -"t + x + ...
2
Newmark obciÄ…Å‚ pozostaÅ‚e wyrazy szeregu, wprowadziÅ‚ współczynniki ² i Å‚ oraz zapisaÅ‚
równania w postaci:
2
("t)
3
& &&t -"t &&&t -"t
xt = xt -"t + "txt -"t + x + ²("t) x
2
2
& & && &&&t -"t
xt = xt -"t + "txt -"t + Å‚ ("t) x
Je\
eli zakładamy, \
e przyspieszenie jest liniowe w przedziale długości kroku, mo\
na
zapisać zale\
ność
&&t &&t
&&&t -"t x - x -1
x =
"t
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 39
Metoda Newmarka
&&&t -"t
x
Po podstawieniu otrzymywane są równania
2
&&
("t) (xt - &&t-1
x )
3
& &&t-"t
xt = xt-"t + "txt-"t + x + ²("t)
2 "t
&&
(xt - &&t-1
x )
2
& & &&
xt = xt-"t + "txt-"t + Å‚ ("t)
"t
Po przekształceniach, otrzymywane są standardowe równania Newmarka
1
ëÅ‚ öÅ‚
2 2
& &&t-"t &&t
xt = xt-"t + "txt-"t + - ² ("t) x + ²("t) x
ìÅ‚ ÷Å‚
2
íÅ‚ Å‚Å‚
& & && &&
xt = xt-"t + (1- Å‚ )"txt-"t + Å‚"txt
Zale\
ności te mo\
na rozwiązać rekurencyjnie dla ka\
dego kroku czasowego, dla ka\
dego
&&t
x
stopnia swobody układu będącego przemieszczeniem. Wyraz wyznacza się z równania
ruchu po podzieleniu jego obu stron przez masÄ™
K Ft
&&t C &
x + xt + xt =
M M M
C K
&&t Ft - xt - xt
&
x =
M M M
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 40
Metoda Newmarka  bezpośrednie rozwiązanie
równań w pojedynczym kroku
Równania Newmarka mo\
na przekształcić do postaci
&&t & &&
x = b1(xt - xt-"t )+ b2xt-"t + b3xt-"t
& & & & &&
xt = b4(xt - xt-"t )+ b5xt-"t + b6xt-"t
gdzie poszczególne współczynniki wynoszą
1 1 1
b1 = , b2 = - , b3 = 1- , b4 = Å‚"tb1, b5 = 1+ Å‚b2"t, b6 = "t(1+ Å‚b3 - Å‚ )
2
²"t 2²
²("t)
Po podstawieniu do równania ruchu układu otrzymuje się zale\
ność
&& &
Mxt + Cxt + Kxt = Ft
& && & &&
(b1M + b4C + K)xt = Ft + M(b1xt-"t - b2xt-"t - b3xt-"t ) + C(b4xt-"t - b5xt-"t - b6xt-"t )
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 41
Stabilność metody Newmarka
W przypadku zerowego tłumienia metoda Newmarka jest warunkowo stabilna je\
eli
1 1 1 Å‚
Å‚ e" , ² d" oraz "t d" - ²
2 2 Émax 2
gdzie Émax jest maksymalnÄ… czÄ™stoÅ›ciÄ… ukÅ‚adu strukturalnego
Metoda Newmarka jest bezwarunkowo stabilna je\
eli spełniony jest warunek
1
2² e" Å‚ e"
2
Jednak w momencie gdy ł jest większa ni\
0,5 wprowadzone zostają błędy związane z
 tłumieniem numerycznym oraz  wydłu\
aniem okresu czasu
Zazwyczaj przyjmuje siÄ™
1 1
Å‚ = , ² =
4 2
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 42
Metoda Newmarka  algorytm obliczeń
Dla układu o jednym stopniu swobody o masie M, sztywności K, tłumieniu C przy wymuszeniu siłą Ft
1. Obliczenia wstępne
1. Przyjąć wartości M, K i C
2. Przyjąć parametry caÅ‚kowania ² oraz Å‚
3. Obliczyć współczynniki (stałe całkowania) b1, b2, b2, b3, b4, b5, b6
4. Wyznaczyć efektywną sztywność
K = K + b1M + b4C
& &&0
5. Przyjąć warunki początkowe: x0, x0, x
t = "t, 2"t, 3"t,...
2. Obliczenia dla ka\
dego kroku czasu:
1. Obliczyć efektywne wymuszenie
& && & &&
Ft = Ft + M(b1xt-"t - b2xt-"t - b3xt-"t ) + C(b4xt-"t - b5xt-"t - b6xt-"t )
2. Wyznaczyć przemieszczenie dla czasu t
Ft
Kxt = Ft Ò! xt =
K
3. Obliczyć przyspieszenie i prędkość dla czasu t
&&t & &&
x = b1(xt - xt-"t )+ b2xt-"t + b3xt-"t
& & & & &&
xt = b4(xt - xt-"t )+ b5xt-"t + b6xt-"t
4. Przejść do kroku 2.1 wraz z t = t + "t
Metody numeryczne i statystyka - Rozwiązywanie równań ró\
niczkowych 43