KATEDRA URZ
DZEC
I SYSTEMÓW AUTOMATYKI
LABORATORIUM
METOD NUMERYCZYCH
W TECHNICE
CAA
KOWANIE NUMERYCZNE
INSTRUKCJA LABORATORYJNA
WERSJA ROBOCZA
OPRACOWAA
:
mgr inż. Tomasz KWAŚNIEWSKI
POLITECHNIKA ÅšWI
TOKRZYSKA
KIELCE 2011
1. Całkowanie numeryczne
Całkowanie numeryczne polega na obliczeniu całki oznaczonej na podstawie funkcji
podcałkowej w pewnych punktach przedziału całkowania [1]. Odpowiednie zależności dające
poszukiwaną wartość przybliżoną całki nazywane są kwadraturami. Funkcję podcałkową
zastępuje się w przedziale [a, b] funkcją interpolującą lub aproksymującą o możliwie prostej
postaci dla której znana jest funkcja pierwotna. Punkty w których obliczane są wartości
funkcji podcałkowej występującej w kwadraturze nazywane są węzłami kwadratury.
Dana jest funkcja f(x) ciągła w przedziale [a, b]:
b
Jeżeli F (x) = f(x) to f (x)dx = F(b) - F(a) F - funkcja pierwotna funkcji f
+"
a
W instrukcji laboratoryjnej rozpatrywane będą trzy kwadratury:
" metoda prostokątów,
" metoda trapezów,
" metoda Simpsona.
2. Metoda prostokątów
W metodzie prostokątów korzystamy z definicji całki oznaczonej Riemanna, w której
wartość całki interpretowana jest jako suma pól obszarów pod wykresem krzywej w zadanym
przedziale całkowania [a, b]. Sumę tę przybliżamy przy pomocy sumy pól odpowiednio
dobranych prostokątów. Sposób postępowania jest następujący:
Rysunek 1. Złożona metoda prostokątów.
Algorytm:
Przedział całkowania [a, b] dzielimy na n równo odległych punktów x1,x2,...,xn. Punkty
wyznaczamy w prosty sposób wg wzoru:
hÅ" 2Å"i -1
( )
xi = a + , dla i = 1,2, & n.
2
b - a
odległość pomiędzy kolejnymi punktami całkowania wyznaczamy z h = ,
n
 obliczamy wartość funkcji w punktach xi,
wartość całki w postaci przybliżonej otrzymujemy sumując pola powierzchni
wszystkich prostokątów.
b
n
(2i
ëÅ‚ -1)h öÅ‚
f (x) dx E" h fi , gdzie fi = f a +
ìÅ‚ ÷Å‚
"
+"
2
íÅ‚ Å‚Å‚
i=1
a
błąd dla metody prostokątów wyraża się wzorem:
1
R = (b - a)Å" h Å" f '(Â ), Â "[a, b]
2
Przykład:
2
Oblicz wartość całki x2 + x +1 dx metodą prostokątów, n = 5.
( )
+"
1
f (x) = x2 + x +1, x" 1, 2 , n=5
[ ]
2 -1
Krok 3. x3 = x2 + h = 1.5
Krok 1. h = = 0.2 ,
5
f3(x3) = x32 + x3 +1 = 4.75
0.2Å" 2Å"1-1
( )
x1 =1+
2
Krok 4. x4 = x3 + h = 1.7
f1(x1) = x12 + x1 +1 = 3.31
f4(x4) = x42 + x4 +1 = 5.59
Krok 2. x2 = x1 + h = 1.3
Krok 5. x5 = x4 + h = 1.9
f2(x2) = x22 + x2 +1 = 3.99
f5(x5) = x52 + x5 +1 = 6.51
2
f (x) dx = h Å" f1 + f2 + f3 + f4 + f5 = 4.83
( )
+"
1
3. Metoda trapezów
Metoda prostokątów przedstawiona w rozdziale powyżej nie jest zbyt dokładna,
dlatego że, pola prostokątów użytych w metodzie z
le odwzorowują pole pod krzywą. Dużo
lepsze odwzorowanie pola pod krzywą otrzymuje się stosując trapezy. Przedział całkowania
[a, b] dzielimy na n+1 równo odległych punktów x0,x1,...,xn. Punkty te wyznaczamy w prosty
sposób wg wzoru:
Rysunek 2. Złożona metoda trapezów.
Algorytm:
b - a
dzielimy przedział całkowania [a, b] na n równych odcinków h = ,
n
wartość funkcji obliczana jest w n+1 węzłach powstałych po podzieleniu przedziału
całkowania,
wartość całki w postaci przybliżonej otrzymujemy sumując pola powierzchni
wszystkich powstałych trapezów.
b
ëÅ‚ fa n fb öÅ‚
ìÅ‚ ÷Å‚
f (x) dx E" hìÅ‚ + fi + , gdzie fa = f1, fb = fn+1
"
+"
÷Å‚
2 2
íÅ‚ i=2 Å‚Å‚
a
błąd dla metody trapezów wyraża się wzorem:
1
R = (b - a)Å" h2 Å" f ''(Â ), Â "[a, b]
12
Przykład:
2
Oblicz wartość całki x2 + x +1 dx metodą trapezów, n = 5.
( )
+"
1
f (x) = x2 + x +1, x" 1, 2 , n=5
[ ]
2 -1
Krok 1. h = = 0.2 ,
5
x1 = a =1
f1(x1) = x12 + x1 +1 = 3
Krok 2. x2 = x1 + h = 1.2
f2(x2) = x22 + x2 +1 = 3.64
Krok 3. x3 = x2 + h = 1.4
f3(x3) = x32 + x3 +1 = 4.36
Krok 4. x4 = x3 + h = 1.6
f4(x4) = x42 + x4 +1 = 5.16
Krok 5. x5 = x4 + h = 1.8
f5(x5) = x52 + x5 +1 = 6.04
Krok 6. x6 = x5 + h = b = 2
f6 (x6) = x62 + x6 +1 = 7
2
ëÅ‚ fa n fb öÅ‚ 3 7
ëÅ‚ öÅ‚
f (x) dx = h + fi + = 0.2Å" + 3.64 + 4.36 + 5.16 + 6.04 + = 4.84
( )
"
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚
+"
2 2 2 2
íÅ‚ Å‚Å‚
íÅ‚ i=2 Å‚Å‚
1
4. Metoda parabol   Simpsona  
Metoda Simpsona jest jedną z dokładniejszych metod przybliżonego całkowania.
W metodzie prostokątów całka oznaczona przybliżana była funkcjami stałymi  liczona była
suma pól prostokątów. W metodzie trapezów całkę przybliżono za pomocą funkcji liniowych
 liczona była suma pól trapezów. W metodzie Simpsona stosujemy jako przybliżenie
parabolę  będziemy obliczać sumy wycinków obszarów pod parabolą. Algorytm jest
następujący:
Rysunek 3. Złożona metoda Simpsona.
Algorytm:
b - a
dzielimy przedział całkowania [a, b] na n równych odcinków; h = ,
n
n
b
2
h
f (x) dx E" f2i-1 + 4 f2i + f2i+1) =
"(
+"
3
i=1
a
h
( f1 + 4( f2 + f4 + ... + fn )+ 2( f3 + f5 + ... + fn-1)+ fn+1)
3
błąd dla metody Simpsona wyraża się wzorem:
1
(4)
R = (b - a)Å" h4 Å" f (Â ), Â "[a, b]
180
Przykład:
2
Oblicz wartość całki x2 + x +1 dx metodą Simpsona, n = 6.
( )
+"
1
f (x) = x2 + x +1, x" 1, 2 , n=6
[ ]
2 -1 1 5
Krok 1. h = = , Krok 5. x5 = x4 + h =
6 6 3
49
x1 = a =1
f5(x5) = x52 + x5 +1 =
9
f1(x1) = x12 + x1 +1 = 3
7 11
Krok 2. x2 = x1 + h = Krok 6. x6 = x5 + h =
6 6
127 223
f2 (x2) = x22 + x2 +1 = f6(x6) = x62 + x6 +1 =
36 36
4
Krok 7. x7 = x6 + h = b = 2
Krok 3. x3 = x2 + h =
3
f7 (x7 ) = x72 + x7 +1 = 7
37
f3(x3) = x32 + x3 +1 =
9
3
Krok 4. x4 = x3 + h =
2
19
f4(x4) = x42 + x4 +1 =
4
n
b
2
h
f (x) dx = f2i-1 + 4Å" f2i + f2i+1 =
( )
"
+"
3
i=1
a
h
f1 + 4Å" f2 + f4 + f6 + 2Å" f3 + f5 + fn+1 =
( ( ) ( ) )
3
1/ 6
( ) ëÅ‚ 127 19 223 37 49 öÅ‚ 29
öÅ‚ öÅ‚
+ 2Å"ëÅ‚ + + 7÷Å‚ = = 4.8 3
( )
ìÅ‚
ìÅ‚3+ 4Å"ëÅ‚ 36 + 4 + 36 ÷Å‚ ìÅ‚ ÷Å‚
3 9 9 6
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚
íÅ‚ Å‚Å‚
5. MATLAB - całkowanie
Program MATLAB oferuje kilka metod numerycznego całkowania funkcji. W
programie mamy możliwość obliczenia całki oznaczonej, jak również wykorzystując
bibliotekę MATLAB a Symbolic Math Toolbox możemy policzyć całkę nieoznaczoną po
uprzednim przekształceniu wyrażeń matematycznych do postaci symbolicznej. Funkcje quad
i quadl wykorzystują metody adaptacyjne. Dzielą przedział całkowania na podprzedziały o
zmiennej długości tak aby najkrótszy przedział wypadał w miejscu największej zmienności
funkcji podcałkowej.
Tabela 1. Całkowanie numeryczne.
Nazwa funkcji Opis funkcji
quad
Całkowanie metodą adaptacyjną Simpsona (3-punktowa reguła Newtona-
Cotes a), dokładna dla wielomianów do 3 stopnia.
quadl
Całkowanie metodą adaptacyjną (4-punktowa reguła Gauss-Lobatto
i 7-punktowa z rozszerzeniem Kronroda) Dokładna dla wielomianów
do 5 i 9 stopnia.
dblquad
Obliczanie całek podwójnych
triplequad
Obliczanie całek potrójnych
trapz,
Całkowanie metodą trapezów  metoda nieadaptacyjna
cumtrapz
Q = quad(fun, a, b, [RelTol, AbsTol], trace)
fun  funkcja podcałkowa,
a, b  granice całkowania,
RelTol, AbsTol  żądana dokładność (względna i bezwzględna),
trace  wypisywanie wartości obliczeń pośrednich.
Przykład:
fun = @(x, a)(a*x.^3 + 2*x.^2 - 7);
tmp = 4
Q = quad(fun, -2, 3, [1.e-5 1.e-3], 1, tmp)
6. Zadania
Zadanie 1.
Napisać funkcje w programie MATLAB rozwiązujące metody całkowania przedstawione
w instrukcji.
Zadanie 2.
Oblicz wartość całki analitycznie oraz przy pomocy funkcji z zadania pierwszego.
6
x2 + 2 dx , n = 10
( )
+"
2
Ä„
Sin(x) dx , n = 8
( )
+"
0
2Å"Ä„
7 Å"Cos(x) - 5x2 + 2 dx , n = 10
( )
+"
0
Narysuj wykres funkcji podcałkowej w przedziale całkowania. Oblicz błąd procentowy dla
każdej z metod całkowania w stosunku do rozwiązania dokładnego oraz błąd samej metody.
Zadanie 3.
Korzystają z metod wbudowanych do całkowania funkcji w programie MATLAB porównać
otrzymane wyniki z rozwiÄ…zaniem otrzymanym w zadaniu poprzednim.
Literatura
[1] J. Povstenko, Wprowadzenie do metod numerycznych, EXIT, Warszawa 2002
[2] Won Y. Yang, Wenwu Cao, Tae S. Chung, John Morris. Applied numerical methods
using MATLAB, John Wiley & Sons, Inc. 2005
KATEDRA URZ
DZEC
I SYSTEMÓW AUTOMATYKI
LABORATORIUM
METOD NUMERYCZYCH
W TECHNICE
RÓWNANIA RÓŻNICZKOWE
INSTRUKCJA LABORATORYJNA
WERSJA ROBOCZA
OPRACOWAA
:
mgr inż. Tomasz KWAŚNIEWSKI
POLITECHNIKA ÅšWI
TOKRZYSKA
KIELCE 2011
1. Równanie różniczkowe
Ogólnie równanie różniczkowe zwyczajne pierwszego rzędu można zapisać:
2
y = f (x, y) (1)
gdzie funkcja f jest daną funkcją dwóch zmiennych. Funkcja y = y(x) jest określona oraz
różniczkowalną w przedziale [a, b] i jest rozwiązaniem równania (1), jeśli
2
y (x) = f (x, y(x)), "x"[a,b] . (2)
Równanie (1) może posiadać wiele rozwiązań. Aby otrzymać jednoznaczność
rozwiązania, do równania (1) trzeba dodać dodatkowe warunki. Mogą to być:
warunki poczÄ…tkowe
y(a) = ya , (3)
lub warunki brzegowe
g(y(a), y(b)) = 0 . (4)
Równanie (1) oraz (3) nazywamy zagadnieniem początkowym lub zagadnieniem Cauchy ego,
a równanie (1) i (4) nazywamy zagadnieniem brzegowym.
Zagadnienie Cauchy ego
Modelowanie matematyczne wielu procesów i zjawisk doprowadza do rozwiązywania
równań różniczkowych. Niektóre rozwiązania równań różniczkowych zwyczajnych mogą być
rozwiązane metodami analitycznymi, ale większość tych równań może być rozwiązana tylko
w sposób numeryczny. W najprostszym przypadku równania różniczkowego zwyczajnego
pierwszego rzędu trzeba znalez
ć funkcję różniczkowalną spełniającą równie (1).
Wiadomo, że rozwiązanie równania (1) nie jest określone jednoznacznie przez samo
równanie. Rozwiązanie ogólne tego równania zależy od dowolnej stałej. Aby otrzymać
rozwiązanie szczególne, potrzebne są dodatkowe warunki. Jednym z najprostszych i
najważniejszych rodzajów zagadnień dla równań różniczkowych zwyczajnych jest
zagadnienie Cauchy ego (zagadnienie początkowe), gdy dodatkowy warunek określony jest w
jednym punkcie. Trzeba więc znalez
ć rozwiązanie równania (1) w przedziale [x0, b]
spełniające warunek początkowy
(2) y(x0) = y0
Twierdzenie: Jeśli funkcja f(x,y) jest ciągła w [x0, b] i spełnia w tym przedziale warunek
Lipschitza
~ ~
(3) f (x, y)- f (x, y) d" L y - y
Gdzie L- stała (stała Lipschitza), to w przedziale (x0, b) istaienje dokładnie jedno rozwiązanie
zagadnienia (1) i (2).
Całkując równanie (1) stronami, otrzymamy równanie całkowe
b
y(x) = y0 + f (x, y)dx
+"
x0
równoważne zagadnieniu (1) i (2).
Przykład metody Eulera
Rozwiązać metodą Eulera podany problem początkowy:
y ' = 2Å" x2 + y , y 0 =1, x " 0,0.5 , h = 0.1. Sprawdz
błąd metody dla rozwiązania
( ) [ ]
dokÅ‚adnego y x = -4 + 5Å"ex - 4Å" x - 2Å" x2 .
( )
Wzór ogólny metody Eulera yi+1 = yi + hÅ" f xi, yi , i = 0,1, 2K
( )
KROK I. y x0 = y0 , czyli y 0 =1, x0 = 0 , i = 0
( ) ( )
y1 = y0 + h Å" f x0, y0 = y0 + h Å" 2Å" x02 + y0 = 1+ 0.1Å" 2Å"02 +1 = 1.1
( )
( ) ( )
KROK II. y1 = 1.1 , x1 = x0 + h = 0.1, i = 1
2
y2 = y1 + hÅ" f x1, y1 =1.1+ 0.1Å" 0.1 +1.1
( ) ( )
(2Å" )=1.212
KROK III. y2 = 1.212 , x2 = x1 + h = 0.2 , i = 2
2
y3 = y2 + hÅ" f x2, y2 =1.212 + 0.1Å" 0.2 +1.212
( ) ( )
(2Å" )=1.3412
KROK IV. y3 = 1.3412 , x3 = x2 + h = 0.3 , i = 3
2
y4 = y3 + hÅ" f x3, y3 = 1.3412 + 0.1Å" 0.3 +1.3412 1.49332
( ) ( )
(2Å" )=
KROK V. y4 = 1.49332 , x4 = x3 + h = 0.4 , i = 4
2
y5 = y4 + hÅ" f x4, y4 =1.49332 + 0.1Å" 0.4 +1.49332
( ) ( )
(2Å" )=1.67465
KROK VI. y5 = 1.67465 , x5 = x4 + h = 0.5 , i = 5
2
y6 = y5 + hÅ" f x5, y5 =1.67465 + 0.1Å" 0.5 +1.67465
( ) ( )
(2Å" )=1.89212
Wyznaczamy błąd procentowy metody Eulera:
y x = -4 + 5Å"ex - 4Å" x - 2Å" x2
( )
0.5
ydok = y x dx = 0.660273
( )
+"
0
ydok + y 0 - y6
( ( ) 0.660273 +1 -1.89212
) ( )
´ = Å"100% = Å"100% =13.96%
ydok + y 0 0.660273+1
( ( ) ( )
)
Rozwiązać metodą Runge-Kutty II rzędu podany problem początkowy:
y ' = 2Å" x2 + y , y 0 =1, x " 0,0.5 , h = 0.1. Sprawdz
błąd metody dla rozwiązania
( ) [ ]
dokÅ‚adnego y x = -4 + 5Å"ex - 4Å" x - 2Å" x2 .
( )
Wzór ogólny metody Runge-Kutty II rzędu
KROK I. y x0 = y0 , czyli y 0 =1, x0 = 0 , i = 0
( ) ( )
g1 = h Å" f x0, y0 = h Å" 2Å" x02 + y0 = 0.1Å" 2Å"02 +1 = 0.1
( )
( ) ( )
2 2
g2 = h Å" f x0 + h, y0 + g1 = h Å" x0 + h + y0 + g1 0.1Å" 0 + 0.1 + 1+ 0.1 0.112
( ) ( ) ( ) ( ) ( )
(2Å" )= (2Å" )=
1 1
y1 = y0 + Å" g1 + g2 = 1+ Å" 0.1+ 0.112 = 1.106
( ) ( )
2 2
KROK II. y1 = 1.106 , x1 = x0 + h = 0.1, i = 1
2
g1 = h Å" f x1, y1 = 0.1Å" 0.1 +1.106 0.1126
( ) ( )
(2Å" )=
2
g2 = h Å" f x1 + h, y1 + g1 = 0.1Å" 0.1+ 0.1 + 1.106 + 0.1126 0.12986
( ) ( ) ( )
(2Å" )=
1 1
y2 = y1 + Å" g1 + g2 = 1.106 + Å" 0.1126 + 0.12986 = 1.22723
( ) ( )
2 2
KROK III. y2 = 1.22723, x2 = x1 + h = 0.2 , i = 2
2
g1 = h Å" f x2, y2 = 0.1Å" 0.2 +1.22723 0.130723
( ) ( )
(2Å" )=
2
g2 = h Å" f x2 + h, y2 + g1 = 0.1Å" 0.2 + 0.1 + 1.22723 + 0.130723 0.153795
( ) ( ) ( )
(2Å" )=
1 1
y3 = y2 + Å" g1 + g2 = 1.22723 + Å" 0.130723 + 0.153795 = 1.36949
( ) ( )
2 2
KROK IV. y3 =1.36949 , x3 = x2 + h = 0.3 , i = 3
2
g1 = h Å" f x3, y3 = 0.1Å" 0.3 +1.36949 0.154949
( ) ( )
(2Å" )=
2
g2 = h Å" f x3 + h, y3 + g1 = 0.1Å" 0.3 + 0.1 + 1.36949 + 0.154949 0.184
( ) ( ) ( )
(2Å" )=
1 1
y4 = y3 + Å" g1 + g2 = 1.36949 + Å" 0.154949 + 0.184 = 1.53896
( ) ( )
2 2
KROK V. y4 = 1.53896 , x4 = x3 + h = 0.4 , i = 4
2
g1 = h Å" f x4, y4 = 0.1Å" 0.4 +1.53896 0.185896
( ) ( )
(2Å" )=
2
g2 = h Å" f x4 + h, y4 + g1 = 0.1Å" 0.4 + 0.1 + 1.53896 + 0.185896 0.222486
( ) ( ) ( )
(2Å" )=
1 1
y5 = y4 + Å" g1 + g2 = 1.53896 + Å" 0.185896 + 0.222486 = 1.74315
( ) ( )
2 2
Wyznaczamy błąd procentowy metody Runge-Kutty II rzędu:
y x = -4 + 5Å"ex - 4Å" x - 2Å" x2
( )
0.5
ydok = y x dx = 0.660273
( )
+"
0
ydok + y 0 - y5
( ( ) 0.660273 +1
) ( )-1.74315
´ = Å"100% = Å"100% = 4.99%
ydok + y 0 0.660273 +1
( ( ) ( )
)
Rozwiązać metodą Runge-Kutty IV rzędu podany problem początkowy:
y ' = 2Å" x2 + y , y 0 =1, x " 0,0.5 , h = 0.1. Sprawdz
błąd metody dla rozwiązania
( ) [ ]
dokÅ‚adnego y x = -4 + 5Å"ex - 4Å" x - 2Å" x2 .
( )
Wzór ogólny metody Runge-Kutty IV rzędu
KROK I. y x0 = y0 , czyli y 0 =1, x0 = 0 , i = 0
( ) ( )
g1 = h Å" f x0, y0 = h Å" 2Å" x02 + y0 = 0.1Å" 2Å"02 +1 = 0.1
( )
( ) ( )
2 2
h g1 ëÅ‚ h g1 0.1 0.1
ëÅ‚ öÅ‚ öÅ‚ ëÅ‚ öÅ‚öÅ‚ ëÅ‚ öÅ‚ ëÅ‚1+ öÅ‚öÅ‚
g2 = h Å" f x0 + , y0 + = h Å"ìÅ‚ 2Å"ëÅ‚ x0 + + y0 + = 0.1Å"ìÅ‚ 2Å"ëÅ‚0 + +
ìÅ‚
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚÷Å‚ ìÅ‚ ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚÷Å‚ = 0.1055
2 2 2 2 2 2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚÷Å‚ íÅ‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚÷Å‚
íÅ‚ Å‚Å‚ Å‚Å‚
2 2
h g2 ëÅ‚ h g2 0.1 0.1055
ëÅ‚ öÅ‚ öÅ‚ ëÅ‚ öÅ‚öÅ‚ ëÅ‚ öÅ‚ ëÅ‚1+ öÅ‚öÅ‚ = 0.105775
g3 = h Å" f x0 + , y0 + = h Å"ìÅ‚ 2Å"ëÅ‚ x0 + + y0 + = 0.1Å"ìÅ‚ 2Å"ëÅ‚ 0 + +
ìÅ‚
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚÷Å‚ ìÅ‚ ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚÷Å‚
2 2 2 2 2 2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚÷Å‚ íÅ‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚÷Å‚
íÅ‚ Å‚Å‚ Å‚Å‚
2 2
g4 = h Å" f x0 + h, y0 + g3 = h Å" x0 + h + y0 + g3 0.1Å" 0 + 0.1 + 1+ 0.105775 0.112578
( ) ( ) ( ) ( ) ( )
(2Å" )= (2Å" )=
1 1
y1 = y0 + Å" g1 + 2Å" g2 + 2Å" g3 + g4 =1+ Å" 0.1+ 2Å"0.1055 + 2Å"0.105775 + 0.112578 = 1.10585
( ) ( )
6 6
KROK II. y1 =1.10585, x1 = x0 + h = 0.1, i = 1
2
g1 = h Å" f x1, y1 = 0.1Å" 0.1 +1.10585 0.112585
( ) ( )
(2Å" )=
2
h g1 ëÅ‚ 0.1 0.112585
ëÅ‚ öÅ‚ öÅ‚ ëÅ‚1.10585 + öÅ‚öÅ‚ = 0.120714
g2 = h Å" f x1 + , y1 + = 0.1Å"ìÅ‚ 2Å"ëÅ‚ 0.1+ +
ìÅ‚
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚÷Å‚
2 2 2 2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚÷Å‚
íÅ‚ Å‚Å‚
2
h g2 ëÅ‚ 0.1 0.120714
ëÅ‚ öÅ‚ öÅ‚ ëÅ‚1.10585 + öÅ‚öÅ‚ = 0.121121
g3 = h Å" f x1 + , y1 + = 0.1Å"ìÅ‚ 2Å"ëÅ‚0.1+ +
ìÅ‚
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚÷Å‚
2 2 2 2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚÷Å‚
íÅ‚ Å‚Å‚
2
g4 = h Å" f x1 + h, y1 + g3 = 0.1Å" 0.1+ 0.1 + 1.10585 + 0.121121 0.130697
( ) ( ) ( )
(2Å" )=
1 1
y2 = y1 + Å" g1 + 2Å" g2 + 2Å" g3 + g4 =1.10585 + Å" 0.112585 + 2Å"0.120714 + 2Å"0.121121+ 0.130697 =1.22701
( ) ( )
6 6
KROK III. y2 = 1.22701, x2 = x1 + h = 0.2 , i = 2
2
g1 = h Å" f x2, y2 = 0.1Å" 0.2 +1.22701 0.130701
( ) ( )
(2Å" )=
2
h g1 ëÅ‚ 0.1 0.130701
ëÅ‚ öÅ‚ öÅ‚ ëÅ‚1.22701+ öÅ‚öÅ‚ = 0.141736
g2 = h Å" f x2 + , y2 + = 0.1Å"ìÅ‚ 2Å"ëÅ‚ 0.2 + +
ìÅ‚
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚÷Å‚
2 2 2 2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚÷Å‚
íÅ‚ Å‚Å‚
2
h g2 ëÅ‚ 0.1 0.141736
ëÅ‚ öÅ‚ öÅ‚ ëÅ‚1.22701+ öÅ‚öÅ‚ = 0.142288
g3 = h Å" f x2 + , y2 + = 0.1Å"ìÅ‚ 2Å"ëÅ‚ 0.2 + +
ìÅ‚
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚÷Å‚
2 2 2 2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚÷Å‚
íÅ‚ Å‚Å‚
2
g4 = h Å" f x2 + h, y2 + g3 = 0.1Å" 0.2 + 0.1 + 1.22701+ 0.142288 0.15493
( ) ( ) ( )
(2Å" )=
1 1
y3 = y2 + Å" g1 + 2Å" g2 + 2Å" g3 + g4 =1.22701+ Å" 0.130701+ 2Å"0.141736 + 2Å"0.142288+ 0.15493 =1.36929
( ) ( )
6 6
KROK IV. y3 =1.36929 , x3 = x2 + h = 0.3 , i = 3
2
g1 = h Å" f x3, y3 = 0.1Å" 0.3 +1.36929 0.154929
( ) ( )
(2Å" )=
2
h g1 ëÅ‚ 0.1 0.154929
ëÅ‚ öÅ‚ öÅ‚ ëÅ‚1.36929 + öÅ‚öÅ‚ = 0.169175
g2 = h Å" f x3 + , y3 + = 0.1Å"ìÅ‚ 2Å"ëÅ‚0.3 + +
ìÅ‚
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚÷Å‚
2 2 2 2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚÷Å‚
íÅ‚ Å‚Å‚
2
h g2 ëÅ‚ 0.1 0.169175
ëÅ‚ öÅ‚ öÅ‚ ëÅ‚1.36929 + öÅ‚öÅ‚ = 0.169888
g3 = h Å" f x3 + , y3 + = 0.1Å"ìÅ‚ 2Å"ëÅ‚0.3 + +
ìÅ‚
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚÷Å‚
2 2 2 2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚÷Å‚
íÅ‚ Å‚Å‚
2
g4 = h Å" f x3 + h, y3 + g3 = 0.1Å" 0.3 + 0.1 + 1.36929 + 0.169888 0.185918
( ) ( ) ( )
(2Å" )=
1 1
y4 = y3 + Å" g1 + 2Å" g2 + 2Å" g3 + g4 =1.36929 + Å" 0.154929 + 2Å"0.169175 + 2Å"0.169888 + 0.185918 =1.53912
( ) ( )
6 6
KROK V. y4 = 1.53912 , x4 = x3 + h = 0.4 , i = 4
2
g1 = h Å" f x4, y4 = 0.1Å" 0.4 +1.53912 0.185912
( ) ( )
(2Å" )=
2
h g1 ëÅ‚ 0.1 0.185912
ëÅ‚ öÅ‚ öÅ‚ ëÅ‚1.53912 + öÅ‚öÅ‚ = 0.203708
g2 = h Å" f x4 + , y4 + = 0.1Å"ìÅ‚ 2Å"ëÅ‚ 0.4 + +
ìÅ‚
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚÷Å‚
2 2 2 2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚÷Å‚
íÅ‚ Å‚Å‚
2
h g2 ëÅ‚ 0.1 0.203708
ëÅ‚ öÅ‚ öÅ‚ ëÅ‚1.53912 + öÅ‚öÅ‚ = 0.204597
g3 = h Å" f x4 + , y4 + = 0.1Å"ìÅ‚ 2Å"ëÅ‚ 0.4 + +
ìÅ‚
ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚ ìÅ‚ ÷Å‚÷Å‚
2 2 2 2
íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚ íÅ‚ Å‚Å‚÷Å‚
íÅ‚ Å‚Å‚
2
g4 = h Å" f x4 + h, y4 + g3 = 0.1Å" 0.4 + 0.1 + 1.53912 + 0.204597 0.224372
( ) ( ) ( )
(2Å" )=
1 1
y5 = y4 + Å" g1 + 2Å" g2 + 2Å" g3 + g4 =1.53912 + Å" 0.185912 + 2Å"0.203708 + 2Å"0.204597 + 0.224372 =1.7436
( ) ( )
6 6
Wyznaczamy błąd procentowy metody Runge-Kutty IV rzędu:
y x = -4 + 5Å"ex - 4Å" x - 2Å" x2
( )
0.5
ydok = y x dx = 0.660273
( )
+"
0
ydok + y 0 - y5
( ( ) 0.660273+1
) ( )-1.7436
´ = Å"100% = Å"100% = 5.01%
ydok + y 0 0.660273+1
( ( ) ( )
)
7. MATLAB  równania różniczkowe
Równania różniczkowe, które mamy możliwość rozwiązywania w programie MATLAB
można podzielić na trzy grupy:
a. równania różniczkowe zwyczajne gdzie szukamy rozwiązania równania
różniczkowego dla zadanego warunku początkowego.
b. równania różniczkowe zwyczajne gdzie szukamy rozwiązania równania
różniczkowego dla zadanych warunków granicznych.
c. równania różniczkowe cząstkowe.
Równania różniczkowe zwyczajne pierwszego rzędu
W rozwiązaniach wielu problemów fizycznych, ekonomicznych wymagana jest znajomość
funkcji y=y(t) przy znajomości funkcji y'=f(t, y) oraz warunków początkowych y(a)=y0,
gdzie a oraz y0 sÄ… liczbami rzeczywistymi a f funkcjÄ™ w postaci jawnej. W takim przypadku
mamy do czynienia z równaniem różniczkowym zwyczajnym pierwszego rzędu (ordinary
differential equations - ODEs).
Tabela 1. Opis funkcji różniczkowych w programie MATLAB.
Nazwa funkcji Opis funkcji
ode23
Rozwiązuje zagadnienie początkowe dla równań różniczkowych
zwyczajnych metodą Runge-Kutta rzędu 2 i 3
ode45
Rozwiązuje zagadnienie początkowe dla równań różniczkowych
zwyczajnych metodą Runge-Kutta rzędu 4 i 5
ode113
Rozwiązuje zagadnienie początkowe dla równań różniczkowych
zwyczajnych metodÄ… Adams-Bashforth-Moulton
ode15s
Metoda oparta na formule numerycznego różniczkowania
ode23s
Metoda oparta na zmodyfikowanej formule Rosenbrocka 2 rzędu
[t, y] = ode45 (fun, tspan, y0,options)
gdzie:
fun -- równanie różniczkowe
tspan -- wartościami czasu, dla którego poszukiwane jest rozwiązanie,
y0 -- jest wektorem, w którym przechowywane są wartości rozwiązania układu w chwili
poczÄ…tkowej,
options -- są ustawiane za pomocą funkcji odeset i pozwalają ingerować w parametry
rozwiązywania równania.
Jeżeli zmienna tspan zawiera więcej niż dwa elementy, to metoda zwraca obliczone
wartości y dla tych elementów. Parametry wyjściowe t i y zawierają wektory wartości do
obliczeń t i obliczone dla nich wartości y .
Przykład:
Rozwiązać równanie różniczkowe zwyczajne w postaci o określonych warunkach
poczÄ…tkowych:
y ' = 2Å" x2 + y , y 0 =1, x " 0,0.5 , h=0.01
( ) [ ]
przy pomocy funkcji ode45.
function test_ode45
t = 0:.01:0.5; % skala czasu
y0 = 1; % warunek poczÄ…tkowy
[t, y] = ode45(@fun, t, y0); % wywołanie funkcji ode45
figure(1)
plot(t, y(:, 1), 'rx')
xlabel('t'); ylabel('y');
grid on
function dy = fun(t, y)
dy = 2*t.^2 + y;
end
end
Rozwiązać równanie różniczkowe drugiego rzędu podane w postaci o określonych warunkach
poczÄ…tkowych:
y'' + 5Å" y' - 4Å" y = sin( 10Å"t ), y(0)= 0, y' (0) = 0 , x "[0, 3], h=0.01
przy pomocy funkcji ode45.
function test_ode45
t = 0:.01:3; % skala czasu
y0 = 0; % warunek poczÄ…tkowy
dy0 = 0; % warunek poczÄ…tkowy
[t, y] = ode45(@fun, t, [y0, dy0]); % wywołanie funkcji ode45
figure(1)
plot(t, y(:, 1), 'rx', t, y(:, 2), 'bo')
xlabel('t'); ylabel('y');
grid on
function dy = fun(t, y)
dy_1 = y(2);
dy_2 = -5*y(2) + 4*y(1) + sin(10*t);
dy = [dy_1; dy_2];
end
end
Rozwiązać metodą Eulera podany problem początkowy:
y ' = 2Å" x2 + y , y 0 =1, x " 0,0.5 , h=0.01. Sprawdz
błąd metody dla rozwiązania
( ) [ ]
dokÅ‚adnego y x = -4 + 5Å"ex - 4Å" x - 2Å" x2
( )
function [t, y] = m_Euler(fun, tspan, y0, N)
if nargin < 4 | N<=0,
N= 100;
end
if nargin < 3,
y0 = 0;
end
h = (tspan(2) - tspan(1))/N; % krok
t = tspan(1) + h*[0:N]'; % przedział czasu
y(1,:) = y0(:)';
for k = 1:N
y(k + 1, :) = y(k, :) +h*feval(fun, t(k), y(k, :));
end
end
- skrypt z wywołaniem funkcji Eulera
function test_Euler
t = [0, 0.5]; % granice czasowe
N = 100; % ilość kroków w przedziale czasowym
y0 = 1; % warunek poczÄ…tkowy
[t, y] = m_Euler(@fun, t, y0, N); % wywołanie funkcji
function dy = fun(t, y)
dy = 2*t.^2 + y;
end
end
8. Zadania
Zadanie 1.
Napisać funkcje w programie MATLAB rozwiązujące równania różniczkowe pierwszego
stopnia.
" Metoda Runge-Kutty drugiego rzędu
1
yi+1 = yi + Å" g1 + g2 , i = 0,1, 2K
( )
2
g1 = hÅ" f xi, yi , g2 = hÅ" f xi + h, yi + g1 , i = 0,1,2K
( ) ( )
" Metoda Runge-Kutty czwartego rzędu
1
yi+1 = yi + Å" g1 + 2Å" g2 + 2Å" g3 + g4 , i = 0,1, 2K
( )
6
h g1
ëÅ‚ öÅ‚,
g1 = h Å" f xi , yi , g2 = h Å" f xi + , yi +
( )
ìÅ‚ ÷Å‚
2 2
íÅ‚ Å‚Å‚
h g2
ëÅ‚ öÅ‚, g4 = h Å" f xi + h, yi + g3
g3 = hÅ" f xi + , yi +
( )
ìÅ‚ ÷Å‚
2 2
íÅ‚ Å‚Å‚
Zadanie 2.
Rozwiązać równanie różniczkowe metodami poznanymi na zajęciach.
" y ' = 5Å" x + 3Å" y , y 0 = 2 , x " 0,1 , h = 0.2 . Sprawdz
błąd metody dla rozwiązania
( ) [ ]
1
dokÅ‚adnego y x = -5 + 23Å"e3Å"x -15Å" x .
( )
( )
9
" y ' = x3 + y , y 0 =1, x " 0,1 , h = 0.2 . Sprawdz
błąd metody dla rozwiązania
( ) [ ]
dokÅ‚adnego y x = -6 + 7Å"ex - 6Å" x - 3Å" x2 - x3 .
( )
" y ' = x2 + 2Å" y , y 0 = 2 , x " 0,1 , h = 0.2 . Sprawdz
błąd metody dla rozwiązania
( ) [ ]
1
dokÅ‚adnego y x = -1+ 9Å"e2Å"x - 2Å" x - 2Å" x2 .
( )
( )
4
Zadanie 3.
Korzystają z metod wbudowanych do rozwiązywania równań różniczkowych w programie
MATLAB porównać otrzymane wyniki z rozwiązaniem otrzymanym w zadaniu poprzednim.