KATEDRA URZĄDZEŃ I SYSTEMÓW AUTOMATYKI
LABORATORIUM
METOD NUMERYCZYCH
W TECHNICE
INTERPOLACJA
INSTRUKCJA LABORATORYJNA DO ĆWICZENIA NR 1.
WERSJA ROBOCZA
OPRACOWAŁ:
mgr inŜ. Tomasz KWAŚNIEWSKI
POLITECHNIKA ŚWIĘTOKRZYSKA
KIELCE 2011
1.
Interpolacja
ZałóŜmy Ŝe dane są wartości funkcji
)
(x
f
na zbiorze punktów
n
x
x
x
...,
,
,
1
0
zwanych
w
ę
złami interpolacji. Zadaniem interpolacji jest wyznaczenie przybli
Ŝ
onych warto
ś
ci funkcji
)
(x
f
zwanej funkcj
ą
interpolowan
ą
w punktach nie b
ę
d
ą
cych w
ę
złami interpolacji.
Przybli
Ŝ
on
ą
warto
ść
funkcji
)
(x
f
obliczamy za pomoc
ą
funkcji
)
( x
F
zwan
ą
funkcj
ą
interpoluj
ą
c
ą
, która w w
ę
złach ma te same warto
ś
ci co funkcja interpolowana. Funkcja
interpoluj
ą
ca jest funkcj
ą
pewnej klasy. Najcz
ęś
ciej b
ę
dzie to wielomian algebraiczny,
wielomian trygonometryczny, funkcja wymierna lub funkcja sklejana. Interpolacj
ę
stosuje si
ę
najcz
ęś
ciej gdy nie znamy analitycznej postaci funkcji
)
(x
f
(jest ona tylko stablicowana) lub
gdy jej posta
ć
analityczna jest zbyt skomplikowana.
2.
Interpolacja wielomianowa, wzory Lagrange’a
Dane s
ą
w
ę
zły interpolacyjne
n
x
x
x
...,
,
,
1
0
, parami ró
Ŝ
ne tzn.
j
i
x
x
j
i
≠
⇔
≠
. Dane
s
ą
warto
ś
ci funkcji interpolowanej w w
ę
złach
n
f
f
f
...,
,
,
1
0
gdzie
)
(
i
i
x
f
f
=
n
i
,...,
2
,
1
=
.
Zadanie interpolacji polega na znalezieniu wielomianu
)
(
x
L
n
stopnia co najwy
Ŝ
ej n, by
warto
ś
ci tego wielomianu i funkcji interpolowanej w w
ę
złach były sobie równe czyli
.
,...,
2
,
1
)
(
n
i
f
x
L
i
i
n
=
=
[1]
Twierdzenie 1
. Zadanie interpolacji wielomianowej posiada jednoznaczne rozwi
ą
zanie, czyli
istnieje tylko jeden wielomian spełniaj
ą
cy warunek [1].
Konstruujemy funkcje pomocnicze:
( )
(
)
(
)
∏
≠
=
−
−
=
n
i
j
j
j
i
j
i
x
x
x
x
x
l
0
funkcje te s
ą
wielomianami stopnia
n takimi,
Ŝ
e
( )
≠
=
=
k
i
dla
k
i
dla
x
l
k
i
0
1
St
ą
d wynika,
Ŝ
e
( )
( ) ( )
( )
(
)
(
)
( ) ( )
i
i
n
i
j
j
j
i
j
n
i
i
i
n
i
i
x
f
x
L
x
x
x
x
x
f
x
l
x
f
x
L
=
⇒
−
−
=
⋅
=
∏
∑
∑
≠
=
=
=
0
0
0
[2]
( ) (
) (
)
(
)
(
) (
)
(
)
(
) (
)
(
)
(
) (
)
(
)
(
) (
)
(
)
(
) (
)
(
)
N
N
N
N
N
N
N
N
N
N
y
x
x
x
x
x
x
x
x
x
x
x
x
y
x
x
x
x
x
x
x
x
x
x
x
x
y
x
x
x
x
x
x
x
x
x
x
x
x
x
L
⋅
−
⋅
⋅
−
⋅
−
−
⋅
⋅
−
⋅
−
+
+
+
⋅
−
⋅
⋅
−
⋅
−
−
⋅
⋅
−
⋅
−
+
⋅
−
⋅
⋅
−
⋅
−
−
⋅
⋅
−
⋅
−
=
−
−
1
1
0
1
1
0
1
1
2
1
0
1
2
0
0
0
2
0
1
0
2
1
K
K
K
K
K
K
K
[3]
Wzór [2] i [3] nosz
ą
nazw
ę
wzoru interpolacyjnego Lagrange’a.
3.
Przykład zadania interpolacji wielomianowej
Stosuj
ą
c metody Lagrange’a zbudowa
ć
wielomian interpolacyjny 4-go stopnia dla
nast
ę
puj
ą
cej tablicy
.
x
0 2 3 5 6
f(x) 1 3 2 7 17
Szukamy wielomianu Lagrange’a stopnia o jeden mniejszego ni
Ŝ
liczba w
ę
złów interpolacji.
W tym przypadku b
ę
dzie to wielomian czwartego stopnia w postaci:
( )
2
3
4
1
2
3
4
5
L x
a
a x
a x
a x
a x
= +
+
+
+
Nast
ę
pnie tworzymy wielomian Lagrange’a korzystaj
ą
c z wzoru [3].
( ) (
) (
) (
) (
)
(
) (
) (
) (
)
(
) (
) (
) (
)
(
) (
) (
) (
)
(
) (
) (
) (
)
(
) (
) (
) (
)
(
) (
) (
) (
)
(
) (
) (
) (
)
(
) (
) (
) (
)
(
) (
) (
) (
)
4
2
3
5
6
0
3
5
6
1
3
0 2
0 3
0 5
0 6
2 0
2 3
2 5
2 6
0
2
5
6
0
2
3
6
2
7
3 0
3 2
3 5
3 6
5 0
5 2
5 3
5 6
0
2
3
5
6 0
6 2
6 3
6 5
x
x
x
x
x
x
x
x
L
x
x
x
x
x
x
x
x
x
x
x
x
x
− ⋅ − ⋅ − ⋅ −
− ⋅ − ⋅ − ⋅ −
=
⋅ +
⋅ +
− ⋅ − ⋅ − ⋅ −
− ⋅ − ⋅ − ⋅ −
− ⋅ − ⋅ − ⋅ −
− ⋅ − ⋅ − ⋅ −
+
⋅ +
⋅ +
− ⋅ − ⋅ − ⋅ −
− ⋅ − ⋅ − ⋅ −
− ⋅ − ⋅ − ⋅ −
+
− ⋅ − ⋅ − ⋅ −
17
⋅
Po
wykonaniu
prostych
przekształce
ń
matematycznych
otrzymujemy
wielomian
interpolacyjny Lagrange’a czwartego stopnia w postaci:
( )
2
3
4
4
47
481
19
1
1
10
180
45
180
L
x
x
x
x
x
= +
⋅ −
⋅ +
⋅ −
⋅
4.
Skrypt wielomianu interpolacyjnego Lagrange’a.
function
[q, L] = lagrange(x, y)
%Input: x = [x0 x1 ... xN]
% y = [y0 y1 ... yN]
%Output: q = współczynniki wielomianu Lagrange'a stopnia N
% L = współczynniki Lagrange'a
N = length(x) - 1;
q = 0;
for
m = 1:N + 1
P = 1;
for
k = 1:N + 1
if
k ~= m,
P = conv(P, [1 -x(k)])/(x(m) - x(k));
end
end
L(m, :) = P;
% wspolczynniki wielomianowe Lagrange'a
q = q + y(m)*P;
% wielomian Lagrange'a
end
5.
Wzór interpolacyjny Newtona
Iloraz róŜnicowy
Funkcja f(x) jest okre
ś
lona za pomoc
ą
tablicy: x
0
, x
1
, …, x
n
s
ą
w
ę
złami interpolacji, a
f(x
0
), f(x
1
), …, f(x
n
) – odpowiadaj
ą
cymi tym w
ę
złom warto
ś
ciami funkcji f(x). Oczywi
ś
cie
x
i
≠
x
j
dla i
≠
j (i, j = 0, 1, …, n) i ponadto ró
Ŝ
nice
∆
x
i
= x
i+1
– x
i
(i = 0, 1, 2, …) nie s
ą
na ogół
stałe.
Wyra
Ŝ
enia:
[
]
( ) ( )
[
]
( ) ( )
[
]
( ) ( )
1
1
1
1
2
1
2
2
1
0
1
0
1
1
0
,
,
,
−
−
−
−
−
=
−
−
=
−
−
=
n
n
n
n
n
n
x
x
x
f
x
f
x
x
f
x
x
x
f
x
f
x
x
f
x
x
x
f
x
f
x
x
f
K
K
K
K
K
K
K
K
K
nazywamy ilorazami ró
Ŝ
nicowymi pierwszego rz
ę
du. Analogicznie definiujemy ilorazy
ró
Ŝ
nicowe drugiego rz
ę
du:
[
]
(
) (
)
[
]
(
) (
)
[
]
(
) (
)
2
1
2
1
1
2
1
3
2
1
3
2
3
2
1
0
2
1
0
2
1
2
1
0
,
,
,
,
,
,
,
,
,
,
,
,
−
−
−
−
−
−
−
−
=
−
−
=
−
−
=
n
n
n
n
n
n
n
n
n
x
x
x
x
f
x
x
f
x
x
x
f
x
x
x
x
f
x
x
f
x
x
x
f
x
x
x
x
f
x
x
f
x
x
x
f
K
K
K
K
K
K
K
K
K
Ogólnie iloraz ró
Ŝ
nicowy rz
ę
du n
tworzymy z ilorazu ró
Ŝ
nicowego n-1 za pomoc
ą
wzoru
rekurencyjnego:
[
]
(
) (
)
i
n
i
n
i
i
i
n
i
i
i
n
i
i
i
x
x
x
x
x
f
x
x
x
f
x
x
x
f
−
−
=
+
−
+
+
+
+
+
+
+
1
1
2
1
1
,
,
,
,
,
,
,
,
,
K
K
K
[4]
dla n = 1, 2, … oraz i = 0, 1, 2, …
Z ilorazów ró
Ŝ
nicowych tworzy si
ę
tablice – oto zasada tworzenia:
Ilorazy ró
Ŝ
nicowe
x
i
f(x
i
)
rz
ę
du 1
rz
ę
du 2
rz
ę
du 3
rz
ę
du 4
x
0
x
1
x
2
x
3
x
4
f(x
0
)
f(x
1
)
f(x
2
)
f(x
3
)
f(x
4
)
f(x
0
, x
1
)
f(x
1
, x
2
)
f(x
2
, x
3
)
f(x
3
, x
4
)
f(x
0
, x
1
, x
2
)
f(x
1
, x
2
, x
3
)
f(x
2
, x
3
, x
4
)
f(x
0
, x
1
, x
2
, x
3
)
f(x
1
, x
2
, x
3
, x
4
)
f(x
0
, x
1
, x
2
, x
3
, x
4
)
Ilorazem róŜnicowym
rz
ę
du k funkcji f opartym na parami ró
Ŝ
nych w
ę
złach x
l
, x
l+1
, …, x
l+k
,
w których okre
ś
lona jest funkcja f nazywamy wyra
Ŝ
enie:
( )
∑
∏
+
=
+
≠
=
+
+
−
=
k
l
l
i
k
l
i
j
l
j
j
i
i
k
l
l
l
x
x
x
f
f
x
x
x
)
(
]
;
,
,
,
[
1
K
[5]
przy czym przez iloraz ró
Ŝ
nicowy rz
ę
du zerowego nazywamy warto
ść
tej funkcji w danym
w
ęź
le: [x
l
; f] = f(x
l
).
Własności ilorazu róŜnicowego
1.
Iloraz ró
Ŝ
nicowy jest funkcj
ą
symetryczn
ą
, czyli
]
;
,
,
,
[
]
;
,
,
,
[
1
1
f
x
x
x
f
x
x
x
k
l
l
l
i
i
i
k
l
l
l
+
+
=
+
+
K
K
gdzie
k
l
l
l
x
x
x
+
+
,
,
,
1
K
jest dowoln
ą
permutacj
ą
liczb l, l+1, …, l+k
2.
Iloraz ró
Ŝ
nicowy jest funkcjonałem liniowym ze wzgl
ę
du na funkcj
ę
f,
tzn. je
ś
li f(x) = g(x) + h(x), to
]
;
,
,
,
[
]
;
,
,
,
[
]
;
,
,
,
[
1
1
1
h
x
x
x
g
x
x
x
f
x
x
x
k
l
l
l
k
l
l
l
k
l
l
l
+
+
+
+
+
+
+
=
K
K
K
3.
Iloraz ró
Ŝ
nicowy dla jednostajnej siatki
ρ
n
(równoodległe w
ę
zły: x
k+1
= x
0
+ (k + 1)h
dla k = 0, 1, …, n–1, gdzie h jest stał
ą
nazywan
ą
długo
ś
ci
ą
kroku):
[
]
( )
n
n
n
h
n
x
f
f
!
;
0
∆
=
ρ
gdzie
n
∆
oznacza ró
Ŝ
nic
ę
progresywn
ą
:
( ) ( )
( )
( ) (
) ( )
( )
( )
(
)
(
) ( )
(
) (
) (
)
(
) ( )
(
)
( )
( )
(
)
0
1
0
0
0
0
0
0
0
0
0
2
0
0
0
0
1
0
0
2
x
f
x
f
x
f
h
x
f
h
x
f
h
x
f
x
f
h
x
f
x
f
x
f
x
f
h
x
f
x
f
x
f
x
f
x
f
k
k
n
−
∆
∆
=
∆
−
+
−
+
−
+
=
−
+
∆
=
∆
∆
=
∆
−
+
=
∆
=
∆
=
∆
Wzór interpolacyjny Newtona
Je
Ŝ
eli teraz podstawimy wzór (5) do wzoru interpolacyjnego Lagrange’a przyjmuj
ą
c,
Ŝ
e l = 0
to otrzymamy nowy wzór nazywany
wzorem interpolacyjnym Newtona
:
( ) ( ) (
) ( ) (
) ( )
(
)
( )
x
x
x
x
f
x
x
x
x
f
x
x
x
f
x
f
x
N
n
n
n
1
1
0
1
2
1
0
0
1
0
0
,
,
,
,
,
,
−
+
+
+
+
=
ω
ω
ω
K
K
[6]
gdzie
( ) (
)(
) (
)
n
n
x
x
x
x
x
x
x
−
−
−
=
K
1
0
ω
[7]
Niech
ρ
n
oznacza dowoln
ą
siatk
ę
bez w
ę
złów wielokrotnych. Wielomianem interpolacyjnym
Newtona N
n
dla funkcji f na siatce
ρ
n
nazywamy wielomian:
( ) ( )
[
]
(
)
[
]
(
)(
)
[
]
(
)(
) (
)
1
1
0
1
0
1
0
2
1
0
0
1
0
0
;
,
,
,
;
,
,
;
,
−
−
−
−
+
+
+
−
−
+
−
+
=
n
n
n
x
x
x
x
x
x
f
x
x
x
x
x
x
x
f
x
x
x
x
x
f
x
x
x
f
x
N
K
K
K
[8]
6.
Skrypt wielomianu interpolacyjnego Newtona.
function
[n, DD] = newton(x, y)
%Input: x = [x0 x1 ... xN]
% y = [y0 y1 ... yN]
%Output: n = współczynniki wielomianu Newtona
N = length(x) - 1;
DD = zeros(N + 1, N + 1);
%w pierwsza kolumn
ę
macierzy kwadratowej N+1 x
N+1
DD(1:N + 1,1) = y';
for
k = 2:N + 1
for
m = 1:N + 2 - k
% tablica ró
Ŝ
nic dzielonych
DD(m, k) = (DD(m + 1, k - 1) - DD(m, k - 1))/(x(m + k - 1) - x(m));
end
end
% okre
ś
lenie współczynników wielomianu interpolacyjnego Newtona
a = DD(1, :);
% równanie (4) instrukcja laboratoryjna
n = a(N + 1);
for
k = N:-1:1
n = [n a(k)] - [0 n*x(k)];
end
5.
Przykład zadania interpolacji metodą Newtona
Korzystaj
ą
c ze wzoru interpolacyjnego Newtona, znale
źć
wielomian interpolacyjny
dla
nast
ę
puj
ą
cej tablicy.
x
0 2 3 5 6
f(x) 1 3 2 7 17
Ilorazy ró
Ŝ
nicowe
x
i
f(x
i
)
rz
ę
du 1 rz
ę
du 2 rz
ę
du 3 rz
ę
du 4
0
2
3
5
6
1
3
2
7
17
1
-1
2,5
10
-2/3
35/30
2,5
11/30
1/3
-1/180
( )
(
) (
)(
)
(
)(
)(
)
(
)(
)(
)(
)
4
3
2
4
4
180
1
45
19
180
481
10
47
1
)
(
5
3
2
0
180
1
3
2
0
30
11
2
0
3
2
0
1
1
x
x
x
x
x
N
x
x
x
x
x
x
x
x
x
x
x
N
−
+
−
+
=
−
−
−
−
−
−
−
−
+
−
−
−
−
+
=
6.
Błąd interpolacji wielomianowej i optymalny dobór węzłów
Jak zaznaczono we wst
ę
pie, za pomoc
ą
wielomianu interpolacyjnego mo
Ŝ
na
wyznaczy
ć
przybli
Ŝ
on
ą
warto
ść
funkcji interpolowanej w punktach nie b
ę
d
ą
cych w
ę
złami.
Je
ś
li funkcja interpolowana
)
(x
f
jest klasy C
n+1
w przedziale domkni
ę
tym
〉
〈 b
a,
oraz
wszystkie w
ę
zły
n
x
x
x
,...,
,
1
0
te
Ŝ
nale
Ŝą
do tego przedziału to bł
ą
d bezwzgl
ę
dny interpolacji
wielomianem Lagrange’a mo
Ŝ
na oszacowa
ć
wzorem
)
(
)!
1
(
)
(
)
(
)
(
1
1
x
n
M
x
L
x
f
x
R
n
n
n
n
+
+
+
≤
−
=
ω
[9]
gdzie
)
(
max
)
1
(
,
1
x
f
M
n
b
a
x
n
+
〉
〈
∈
+
=
a
)
(
1
x
n
+
ω
jest wielomianem czynnikowym okre
ś
lonym wzorem [7].
Zauwa
Ŝ
my,
Ŝ
e bł
ą
d interpolacji zale
Ŝ
y nie tylko od (n+1)-szej pochodnej funkcji
interpolowanej ale równie
Ŝ
od wielomianu czynnikowego
)
(
1
x
n
+
ω
, który z kolei zale
Ŝ
y od
doboru w
ę
złów interpolacji
n
x
x
x
,...,
,
1
0
. Problem optymalnego doboru w
ę
złów
interpolacyjnych rozwi
ą
zał Czebyszew. Wykazał,
Ŝ
e warto
ś
ci w
ę
złów w przedziale
〉
〈 b
a,
,
optymalnie dobranych s
ą
okre
ś
lone
n
i
n
i
a
b
b
a
x
i
,...
1
,
0
)
2
2
1
2
cos(
)
(
2
1
)
(
2
1
=
+
+
−
+
+
=
π
.
[10]
Wtedy najmniejsze oszacowanie bł
ę
du interpolacji wynosi
1
2
1
1
2
)
(
)!
1
(
)
(
)
(
)
(
+
+
+
−
+
≤
−
=
n
n
n
n
n
a
b
n
M
x
L
x
f
x
R
[11]
Optymalnie dobrane w
ę
zły wcale nie s
ą
równo odległe lecz zag
ę
szczaj
ą
si
ę
przy ko
ń
cach
przedziału.
7.
Iteracyjna metoda Aitkena
Wielomian w postaci wzoru Lagrange’a jest niewygodny zarówno do wyznaczania
jego warto
ś
ci w dowolnym punkcie (stosuje si
ę
wzór Aitkena) jak i jego całkowania b
ą
d
ź
ró
Ŝ
niczkowania. Cz
ęś
ciej wielomian interpolacyjny okre
ś
la si
ę
w postaci wzoru Newtona
przy czym obydwa wzory s
ą
sobie równowa
Ŝ
ne poniewa
Ŝ
, zgodnie z twierdzeniem 1, istnieje
tylko jeden wielomian interpolacyjny dla w
ę
złów
n
x
x
x
,...,
,
1
0
.
Istnieje metoda obliczania warto
ś
ci wielomianu Lagrange'a w zadanym punkcie, bez
obliczania samego wielomianu interpolacyjnego. Słu
Ŝ
y do tego iteracyjna metoda Aitkena.
Oznaczmy przez
,
i j
W wielomian który w w
ę
złach
(
)
,
i
j
x x i
j
≠
przyjmuje warto
ś
ci
( )
( )
,
j
i
f x
f x :
,
i
i
j
j
i j
j
i
y
x
x
y
x
x
W
x
x
−
−
=
−
Co mo
Ŝ
na uogólni
ć
jako:
( )
1,2,
,
1,
1,2,
,
1,
1,2,
, ,
k
k
k
k
m
m
k m
m
k
W
x
x
W
x
x
W
x
x
x
−
−
−
−
=
−
K
K
K
Aby obliczy
ć
warto
ść
wielomianu interpolacyjnego opartego na n w
ę
złach w dowolnym
punkcie a ró
Ŝ
nym od w
ę
złów, nale
Ŝ
y obliczy
ć
warto
ść
1,2,
,n
W
K
. Wszystkie wyniki
niezb
ę
dnych oblicze
ń
wygodnie jest umie
ś
ci
ć
w macierzy trójk
ą
tnej wraz z w
ę
złami oraz ich
warto
ś
ciami (schemat taki nazywamy schematem Aitkena). Rozwi
ą
zanie takie jest dogodne
zarówno podczas rachunków r
ę
cznych, jak i maszynowych, gdy
Ŝ
podczas obliczania ka
Ŝ
dej
warto
ś
ci zawsze korzystamy z warto
ś
ci poło
Ŝ
onych na lewo w tym samym rz
ę
dzie i
powy
Ŝ
szych.
1
1
2
2
1,2
3
3
1,3
1,2,3
1,
1,2,
1,2, ,
n
n
n
n
n
x
y
x
y
W
x
y
W
W
x
y
W
W
W
K
K
K
K
K
Wyznaczenie
f(4) metod
ą
Aitkena :
1,2
1
0 4
3
2 4
5
2 0
W
−
−
=
=
−
,
1,3
1
0 4
2
3 4
7
3 0
3
W
−
−
=
=
−
,
1,4
1
0 4
7
5 4
29
5 0
5
W
−
−
=
=
−
,
1,5
1
0 4
17
6 4
70
6 0
6
W
−
−
=
=
−
1,2,3
5
2 4
7
3 4
1
3
3 2
3
W
−
−
=
= −
−
,
1,2,4
5
2 4
29
5 4
83
5
5 2
15
W
−
−
=
=
−
,
1,2,5
5
2 4
70
6 4
25
6
6 2
3
W
−
−
=
=
−
1,2,3,4
1
3 4
3
83
5 4
13
15
5 3
5
W
−
−
−
=
=
−
,
1,2,3,5
1
3 4
3
25
6 4
23
3
6 3
9
W
−
−
−
=
=
−
1,2,3,4,5
13
5 4
5
23
6 4
119
9
6 5
45
W
−
−
=
=
−
.
0
1
2
3
5
7
1
3
2
3
3
29
83
13
5
7
5
15
5
70
25
23
119
6 17
6
3
9
45
−
St
ą
d
1,2,3,4,5
119
45
W
=
, jest warto
ś
ci
ą
funkcji Lagrange’a w punkcie
4
x
=
.
8.
MATLAB - rodzaje interpolacji
W programie MATLAB mamy do dyspozycji kilka metod interpolacji:
- wielomian pierwszego stopnia,
- wielomian trzeciego stopnia,
- metoda najbli
Ŝ
szych s
ą
siadów,
- funkcje sklejane.
Interpolacje stosuje si
ę
do zag
ę
szczania tabel. Zabieg ten pozwala na dokładniejsze
przybli
Ŝ
enie danych zawartych w tabelach. Zadanie interpolacji mo
Ŝ
emy podzieli
ć
na:
Interpolacja nieparametryczna
- algorytm najbli
Ŝ
szy s
ą
siad -
nearest
Interpolacja parametryczna
- wielomianowa –
linear
,
cubic
,
pchip
- funkcjami sklejanymi -
spline
Program MATLAB, w podstawowej bibliotece posiada funkcj
ę
o nazwie
interp1
.
Funkcja umo
Ŝ
liwia rozwi
ą
zanie zadania interpolacji funkcji jednej zmiennej y = f(x)
w punktach x
i
nie b
ę
d
ą
cych w
ę
złami interpolacyjnymi.
yi = interp1(x, y, xi, 'metoda')
gdzie:
x
,
y
- wektory współrz
ę
dnych w
ę
złów interpolacji,
xi
- wektor punktów na osi X dla których b
ę
d
ą
obliczane interpolowane warto
ś
ci
yi
.
9.
Zadania
Zadanie 1.
Znale
źć
wielomian interpolacyjny Lagrange’a dla funkcji okre
ś
lonej tablic
ą
warto
ś
ci oraz
wyznaczy
ć
f(7) wykorzystuj
ą
c metod
ę
AITKENA.
x
-2 1 2
4 9
f(x) 3
1 -3 8 15
wyznaczy
ć
f(7) wykorzystuj
ą
c metod
ę
AITKENA
x
-6 -2 2
4
8
f(x) 2
5
10 12 13
wyznaczy
ć
f(-4) wykorzystuj
ą
c metod
ę
AITKENA
x
-5 -3 2 7 12
f(x) -2 -1 0 8 10
wyznaczy
ć
f(5) wykorzystuj
ą
c metod
ę
AITKENA
Zadanie 2.
Dokona
ć
interpolacji funkcji:
a) -
( )
( )
x
x
f
⋅
=
π
cos
b) -
( )
( )
x
x
x
f
⋅
=
π
cos
3
c) -
( )
( )
2
sin
x
x
f
⋅
=
π
w przedziale
6
,
2
−
∈
x
z krokiem 0,5 funkcj
ą
wbudowan
ą
w Matlab’a interp1 wszystkimi
metodami. Narysowa
ć
wykres danej funkcji i funkcji interpoluj
ą
cej na jednym układzie
współrz
ę
dnych oraz wykres bł
ę
du interpolacji.
Zadanie 3.
Oceni
ć
z jak
ą
dokładno
ś
ci
ą
mo
Ŝ
na obliczy
ć
warto
ść
ln 100,5 przy u
Ŝ
yciu wzoru
interpolacyjnego Lagrange’a, je
Ŝ
eli dane s
ą
warto
ś
ci: ln 100, ln 101, ln 102 i ln 103.
Zadanie 4.
Z jak
ą
dokładno
ś
ci
ą
mo
Ŝ
na obliczy
ć
warto
ść
36
sin
π
stosuj
ą
c wzór interpolacyjny Lagrange’a,
je
Ŝ
eli dane s
ą
warto
ś
ci
3
sin
,
4
sin
,
6
sin
,
0
sin
π
π
π
.
Zadanie 5.
Wykorzystuj
ą
c ilorazy ró
Ŝ
nicowe podaj wzór ogólny wielomianu interpolacyjnego Newtona
dla funkcji okre
ś
lonej tablic
ą
w zadaniu 1.
Zadanie 6.
Obliczy
ć
warto
ść
119
za pomoc
ą
wzoru interpolacyjnego Lagrange,a dla funkcji
x
y
=
i
w
ę
złów interpolacji:
.
144
,
121
,
100
2
1
0
=
=
=
x
x
x
Oszacowa
ć
bł
ą
d.
Program ćwiczenia
1.
Napisa
ć
program, który oblicza warto
ś
ci wielomianu Lagrange’a dla dowolnych
punktów x, w
ę
złów równoodległych i zadanej przez prowadz
ą
cego funkcji
interpolowanej.
2.
Napisa
ć
program, który oblicza warto
ś
ci wielomianu Newtona dla dowolnych
punktów x, w
ę
złów równoodległych i zadanej przez prowadz
ą
cego funkcji
interpolowanej.
3.
Napisa
ć
program, który oblicza warto
ś
ci wielomianu wg schematu Aitkena dla
dowolnych punktów x, w
ę
złów równoodległych i zadanej przez prowadz
ą
cego
funkcji interpolowanej.
4.
Sprawozdanie powinno zawiera
ć
:
a)
kod programu i wyniki przeprowadzonych oblicze
ń
,
b)
opis przeprowadzonych bada
ń
,
c)
analiz
ę
uzyskanych wyników,
d)
wnioski.
7. Pytania kontrolne
1.
Co to jest interpolacja i po co si
ę
j
ą
stosuje?
2.
Sformułuj zadanie interpolacji wielomianowej.
3.
Podaj ogóln
ą
posta
ć
wzoru Lagrange’a. Znale
źć
wielomian interpolacyjny Lagrange’a dla
funkcji okre
ś
lonej przez prowadz
ą
cego (zastosuj wzór Lagrange’a).
4.
Podaj wzór Newtona. Znale
źć
wielomian interpolacyjny Newtona dla funkcji okre
ś
lonej
przez prowadz
ą
cego (zastosuj wzór Newtona).
5.
Co to jest iloraz ró
Ŝ
nicowy stopnia n-tego. Jak wygl
ą
da zale
Ŝ
no
ść
rekurencyjna ilorazu
ró
Ŝ
nicowego?
6.
Podaj wzór na bł
ą
d interpolacji wielomianowej. Jak wygl
ą
da optymalny dobór w
ę
złów
interpolacji?
7.
Scharakteryzuj algorytm Aitkena, po co si
ę
go stosuje?
KATEDRA URZĄDZEŃ I SYSTEMÓW AUTOMATYKI
LABORATORIUM
METOD NUMERYCZYCH
W TECHNICE
APROKSYMACJA
INSTRUKCJA LABORATORYJNA DO
Ć
WICZENIA NR 2.
WERSJA ROBOCZA
OPRACOWAŁ:
mgr in
Ŝ
. Tomasz KWA
Ś
NIEWSKI
POLITECHNIKA ŚWIĘTOKRZYSKA
KIELCE 2011
5.
Aproksymacja
Aproksymacja jest to przybli
Ŝ
anie funkcji
)
(x
f
zwanej funkcj
ą
aproksymowan
ą
inn
ą
funkcj
ą
)
( x
Q
zwan
ą
funkcj
ą
aproksymuj
ą
c
ą
. Aproksymacja bardzo cz
ę
sto wyst
ę
puje w
dwóch przypadkach: gdy funkcja aproksymowana jest przedstawiona w postaci tablicy
warto
ś
ci i poszukujemy dla niej odpowiedniej funkcji ci
ą
głej lub gdy funkcj
ę
o dosy
ć
skomplikowanym zapisie analitycznym chcemy przedstawi
ć
w „prostszej” postaci.
Dokonuj
ą
c aproksymacji funkcji
)
(x
f
musimy rozwi
ą
za
ć
dwa wa
Ŝ
ne problemy [1].
•
Dobór odpowiedniej funkcji aproksymuj
ą
cej
)
( x
Q
. Najcz
ęś
ciej b
ę
dzie to tzw.
wielomian uogólniony b
ę
d
ą
cy kombinacj
ą
liniow
ą
funkcji bazowych
)
( x
q
j
∑
=
=
m
j
j
j
m
x
q
a
x
Q
0
)
(
)
(
[1]
Jako funkcje bazowe stosowane s
ą
:
o
jednomiany
m
x
x
x
,...,
,
,
1
2
,
o
funkcje trygonometryczne
mx
mx
x
x
x
x
cos
,
sin
,...,
2
cos
,
2
sin
,
cos
,
sin
,
1
,
o
wielomiany ortogonalne.
Przyj
ę
cie odpowiednich funkcji bazowych powoduje,
Ŝ
e aby wyznaczy
ć
funkcj
ę
aproksymuj
ą
c
ą
nale
Ŝ
y wyznaczy
ć
warto
ś
ci współczynników
m
a
a
a
,...,
,
1
0
.
•
Okre
ś
lenie dokładno
ś
ci dokonanej aproksymacji. Aproksymacja funkcji powoduje
powstanie bł
ę
dów i sposób ich oszacowania wpływa na wybór metody aproksymacji.
Je
ś
li bł
ą
d b
ę
dzie mierzony na dyskretnym zbiorze punktów
n
x
x
x
,...,
,
1
0
to jest to
aproksymacja punktowa,
a je
ś
li b
ę
dzie mierzony w przedziale
〉
〈 b
a
,
to jest to
aproksymacja integralna
lub przedziałowa.
Najcz
ęś
ciej stosowanymi miarami bł
ę
dów aproksymacji s
ą
:
o
dla aproksymacji
ś
redniokwadratowej punktowej
∑
=
−
=
n
i
i
i
x
Q
x
f
S
0
2
)}
(
)
(
{
o
dla aproksymacji
ś
redniokwadratowej integralnej lub przedziałowej
∫
−
=
b
a
dx
x
Q
x
f
S
2
)}
(
)
(
{
o
dla aproksymacji jednostajnej
)
(
)
(
sup
,
x
Q
x
f
S
b
a
x
−
=
〉
〈
∈
.
We wszystkich tych przypadkach zadanie aproksymacji sprowadza si
ę
do takiego
optymalnego doboru funkcji aproksymuj
ą
cej (dla wielomianów uogólnionych za
ś
do
optymalnego doboru współczynników
m
a
a
a
,...,
,
1
0
) aby zdefiniowane wy
Ŝ
ej bł
ę
dy były
minimalne.
6.
Aproksymacja średniokwadratowa
Dane s
ą
punkty
n
x
x
x
,...,
,
1
0
parami ró
Ŝ
ne czyli
j
i
x
x
j
i
≠
⇔
≠
oraz dane s
ą
warto
ś
ci
funkcji aproksymowanej w tych punktach
n
o
f
f
f
,...,
,
1
gdzie
)
(
i
i
x
f
f
=
n
i
,...,
1
,
0
=
.
Zadaniem aproksymacji jest znale
źć
warto
ś
ci współczynników
m
a
a
a
,...,
,
1
0
wielomianu
)
( x
Q
m
stopnia m-tego o postaci
∑
=
=
m
j
j
j
m
x
a
x
Q
0
)
(
[2]
aby bł
ą
d
ś
redniokwadratowy był najmniejszy czyli
∑
∑
=
=
−
=
n
i
m
j
j
i
j
i
a
a
a
x
a
f
S
m
0
2
0
,...,
,
)
(
min
1
0
.
[3]
Zauwa
Ŝ
my,
Ŝ
e bł
ą
d
ś
redniokwadratowy jest funkcj
ą
współczynników
m
a
a
a
,...,
,
1
0
a wi
ę
c z
warunku koniecznego na ekstremum funkcji wielu zmiennych
0
=
∂
∂
j
a
S
m
j
,...,
1
,
0
=
[4]
otrzymamy układ
1
+
m
równa
ń
liniowych o
1
+
m
niewiadomych współczynnikach
m
a
a
a
,...,
,
1
0
, który mo
Ŝ
na zapisa
ć
w postaci sumy
k
m
j
j
j
k
r
a
g
=
∑
=
0
m
k
,...,
1
,
0
=
[5]
gdzie
∑
=
+
=
n
i
j
k
i
j
k
x
g
0
∑
=
=
n
i
k
i
i
k
x
f
r
0
.
[6]
Zadanie aproksymacji
ś
redniokwadratowej punktowej sprowadza si
ę
do rozwi
ą
zania
1
+
m
równa
ń
o
1
+
m
niewiadomych. Rozwa
Ŝ
my jeszcze problem doboru liczby punktów
1
+
n
i
stopnia wielomianu
m
. Je
ś
li n=m to mamy funkcj
ę
aproksymowan
ą
okre
ś
lon
ą
w
1
+
n
punktach i poszukujemy wielomianu aproksymuj
ą
cego stopnia
n
. Jest to wi
ę
c interpolacja i
wtedy
)
(
)
(
i
n
i
x
Q
x
f
=
czyli bł
ą
d
ś
redniokwadratowy
0
=
S
.
Natomiast je
ś
li
m
n
>
, to oznacza to,
Ŝ
e mamy wi
ę
cej punktów ni
Ŝ
wynosi stopie
ń
wielomianu aproksymuj
ą
cego. W tej sytuacji wielomian aproksymuj
ą
cy nie b
ę
dzie dokładnie
odtwarzał wszystkich punktów funkcji aproksymowanej.
7.
Przykład zadania aproksymacji
Znale
źć
wielomian aproksymacyjny zbiór punktów z tabeli poni
Ŝ
ej dla nast
ę
puj
ą
cych funkcji
bazowych
0
1
a
=
,
( )
1
sin
a
x
=
.
i
0
1
2
3
x
6
π
−
0
6
π
2
π
f(x)
1
-8 -3
5
Funkcje bazowe :
( )
0
1
1
sin
a
a
x
=
=
zatem wielomian aproksymuj
ą
cy zbiór punktów ma posta
ć
:
( )
( )
0
0
1
1
0
1
sin
Q x
b a
b a
b
b
x
= ⋅ + ⋅ = + ⋅
Wykorzystuj
ą
c wzór [3] otrzymujemy :
( )
( )
(
)
(
)
(
)
2
3
2
0
1
0
2
2
2
0
0
1
0
1
1
1
2
1
8
3
5
2
i
i
i
i
i
S
Q x
f x
b
b
b
b
b
b
b
=
=
−
=
− ⋅ −
+
+
+
+
+ ⋅ +
+
+ −
∑
Warunek konieczny
0
=
∂
∂
j
b
S
0
1
0
0
1
1
8
2
10
0
2
3
6
0
S
b
b
b
S
b
b
b
∂ = ⋅ + ⋅ + =
∂
∂ = ⋅ + ⋅ − =
∂
Wykonuj
ą
c proste przekształcenia matematyczne otrzymujemy nast
ę
puj
ą
cy wielomian
aproksymuj
ą
cy zbiór punktów podanych w tabeli.
( )
( )
2.1 3.4 sin
Q x
x
= −
+
⋅
10.
MATLAB - aproksymacja
Je
ś
li liczba zadanych punktów jest wi
ę
ksza od stopnia wielomiany
n +1
to cz
ę
sto nie
mo
Ŝ
na przeprowadzi
ć
krzywej przez wszystkie zadane punkty. Mamy wtedy do czynienia z
zadaniem aproksymacji. Rozwi
ą
zanie mo
Ŝ
e by
ć
wielomian, który najlepiej przybli
Ŝ
a zadan
ą
krzyw
ą
w sensie minimum kwadratu bł
ę
du
.
W programie MATLAB zadanie aproksymacji
rozwi
ą
zuje funkcja
polyfit
, która wyznacza współczynniki poszukiwanego wielomianu
aproksymuj
ą
cego:
p = polyfit (x, y, n)
Warto
ś
ci wielomianu w dowolnych punktach zapisanych w wektorze
x1
, oblicza funkcja
polyval
.
y_p = polyval(p, x1)
Parametry wej
ś
ciowe obu funkcji s
ą
nast
ę
puj
ą
ce:
x, y
– odpowiednio wektory warto
ś
ci zmiennej niezale
Ŝ
nej i zale
Ŝ
nej,
n
– stopie
ń
wielomianu aproksymuj
ą
cego,
p
– wektor współczynników poszukiwanego wielomianu,
y_p
– wektor warto
ś
ci wielomianu w dowolnych punktach x1.
11.
MATLAB – Basic Fitting
W programie MATLAB zagadnienie interpolacji lub aproksymacji mo
Ŝ
na rozwi
ą
za
ć
miedzy
innymi wykorzystuj
ą
c okno interfejsu Basic Fitting. Dost
ę
p do tego interfejsu uzyskuje si
ę
poprzez okno graficzne. Na pocz
ą
tku nale
Ŝ
y narysowa
ć
wykres funkcji (interpolowanej lub
aproksymowanej) przy u
Ŝ
yciu dost
ę
pnych funkcji w programie MATLAB [2]. W oknie
graficznym w którym otrzymamy wykres funkcji, nale
Ŝ
y wybra
ć
z menu opcj
ę
Tools -> Basic
Fitting. Uruchomione zostaje okno o nazwie Basic Fitting, które daje nam du
Ŝ
e mo
Ŝ
liwo
ś
ci
wykonywania ró
Ŝ
nego rodzaju interpolacji i aproksymacji zadanych warto
ś
ci funkcji.
Rysunek 1. Interpolacja i aproksymacja z wykorzystaniem interfejsu Basic Fitting.
7.
Zadania
Zadanie 1.
•
Znale
źć
wielomian aproksymacyjny dla funkcji okre
ś
lonej tablic
ą
warto
ś
ci dla funkcji
bazowych w postaci:
0
1
a
=
,
( )
1
cos
a
x
=
.
i
0
1
2
3
x
3
π
−
0
3
π
π
f(x)
3
-1 -3
9
•
Znale
źć
wielomian aproksymacyjny dla funkcji okre
ś
lonej tablic
ą
warto
ś
ci dla funkcji
bazowych w postaci:
0
1
a
=
,
( )
1
sin
a
x
=
.
i
0
1
2
3
x
12
π
−
0
4
π
2
3
π
f(x)
-2
-1
0
2
•
Znale
źć
wielomian aproksymacyjny dla funkcji okre
ś
lonej tablic
ą
warto
ś
ci dla funkcji
bazowych w postaci:
0
1
a
=
,
1
a
x
=
,
2
2
a
x
=
.
i
0
1 2 3
x
-5 -4 5 9
f(x) 12 5 1 3
Zadanie 2.
Znale
źć
najlepszy wielomian aproksymacyjny dla funkcji okre
ś
lonej tablic
ą
warto
ś
ci stosuj
ą
c
okno interfejsu – Basic Fitting.
x
0 1 2 3 4 5 6 7 8 9
f(x) 3 3 3 3 3 6 6 6 6 6
Zadanie 3.
Dokona
ć
aproksymacji
ś
redniokwadratowej funkcji:
•
( )
2
2
+
=
x
x
x
f
wielomianem 2-go stopnia w przedziale
2
,
2
−
∈
x
z krokiem
∆
= 0.01.
•
( )
x
x
x
x
f
−
+
−
=
4
5
2
3
wielomianem 3-go stopnia w przedziale
2
,
2
−
∈
x
z krokiem
∆
= 0.01.
•
( )
2
5
4
2
3
+
−
=
x
x
x
x
f
wielomianem 4-go stopnia w przedziale
2
,
2
−
∈
x
z krokiem
∆
= 0.01.
Narysowa
ć
wykres danej funkcji i funkcji aproksymuj
ą
cej w jednym układzie współrz
ę
dnych,
za
ś
wykres bł
ę
du aproksymacji w drugim. Wyznaczy
ć
maksymaln
ą
warto
ść
bezwzgl
ę
dnego
bł
ę
du aproksymacji w zadanym przedziale.
4.
Pytania kontrolne
1.
Co to jest aproksymacja i po co si
ę
j
ą
stosuje?
2.
Sformułuj problem aproksymacji
ś
redniokwadratowej, punktowej, wielomianowej.
3.
Podaj rozwi
ą
zanie problemu aproksymacji
ś
redniokwadratowej, punktowej,
wielomianowej.
4.
Co to jest wygładzanie funkcji?
Literatura
[1]
A. Jastriebow, M. Wci
ś
lik, Wst
ę
p do metod numerycznych, Politechnika
Ś
wi
ę
tokrzyska, Kielce 2000.
[2]
B. Mrozek, Z. Mrozek, MATLAB i Simulink – poradnik u
Ŝ
ytkownika, Helion 2004
KATEDRA URZĄDZEŃ I SYSTEMÓW AUTOMATYKI
LABORATORIUM
METOD NUMERYCZYCH
W TECHNICE
CAŁKOWANIE NUMERYCZNE
INSTRUKCJA LABORATORYJNA DO
Ć
WICZENIA NR 3.
WERSJA ROBOCZA
OPRACOWAŁ:
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]:
Je
Ŝ
eli
F’(x) = f(x) to
( )
( )
( )
b
a
f x dx
F b
F a
=
−
∫
F - funkcja pierwotna funkcji f
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 x
1
,x
2
,...,x
n
. Punkty
wyznaczamy w prosty sposób wg wzoru:
(
)
2
1
2
i
h
i
x
a
⋅ ⋅ −
= +
, dla i = 1,2, … n.
odległo
ść
pomi
ę
dzy kolejnymi punktami całkowania wyznaczamy z
n
a
b
h
−
=
,
obliczamy warto
ść
funkcji w punktach x
i
,
warto
ść
całki w postaci przybli
Ŝ
onej otrzymujemy sumuj
ą
c pola powierzchni
wszystkich prostok
ą
tów.
−
+
=
≅
∫
∑
=
2
)
1
2
(
,
)
(
1
h
i
a
f
f
gdzie
f
h
dx
x
f
i
b
a
n
i
i
bł
ą
d dla metody prostok
ą
tów wyra
Ŝ
a si
ę
wzorem:
(
)
( )
[ ]
b
a
f
h
a
b
R
,
,
'
2
1
∈
⋅
⋅
−
=
ς
ς
Przykład:
Oblicz warto
ść
całki
(
)
2
2
1
1
x
x
dx
+ +
∫
metod
ą
prostok
ą
tów, n = 5.
2
( )
1
f x
x
x
=
+ +
,
[ ]
1, 2
x
∈
,
n=5
Krok 1.
2 1
0.2
5
h
−
=
=
,
(
)
1
0.2 2 1 1
1
2
x
⋅ ⋅ −
= +
2
1
1
1
1
( )
1 3.31
f x
x
x
=
+ + =
Krok 2.
2
1
1.3
x
x
h
= + =
2
2
2
2
2
(
)
1 3.99
f x
x
x
=
+ + =
Krok 3.
3
2
1.5
x
x
h
= + =
2
3
3
3
3
( )
1
4.75
f x
x
x
=
+ + =
Krok 4.
4
3
1.7
x
x
h
= + =
2
4
4
4
4
(
)
1 5.59
f x
x
x
=
+ + =
Krok 5.
5
4
1.9
x
x
h
= + =
2
5
5
5
5
( )
1
6.51
f x
x
x
=
+ + =
(
)
2
1
2
3
4
5
1
( )
4.83
f x dx
h
f
f
f
f
f
= ⋅
+ + + +
=
∫
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
ź
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 x
0
,x
1
,...,x
n
. Punkty te wyznaczamy w prosty
sposób wg wzoru:
Rysunek 2. ZłoŜona metoda trapezów.
Algorytm:
dzielimy przedział całkowania [a, b] na n równych odcinków
n
a
b
h
−
=
,
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.
1
1
2
,
,
2
2
)
(
+
=
=
=
+
+
≅
∫
∑
n
b
a
b
a
b
n
i
i
a
f
f
f
f
gdzie
f
f
f
h
dx
x
f
bł
ą
d dla metody trapezów wyra
Ŝ
a si
ę
wzorem:
(
)
( )
[ ]
b
a
f
h
a
b
R
,
,
''
12
1
2
∈
⋅
⋅
−
=
ς
ς
Przykład:
Oblicz warto
ść
całki
(
)
2
2
1
1
x
x
dx
+ +
∫
metod
ą
trapezów, n = 5.
2
( )
1
f x
x
x
=
+ +
,
[ ]
1, 2
x
∈
,
n=5
Krok 1.
2 1
0.2
5
h
−
=
=
,
1
1
x
a
= =
2
1
1
1
1
( )
1
3
f x
x
x
=
+ + =
Krok 2.
2
1
1.2
x
x
h
= + =
2
2
2
2
2
(
)
1 3.64
f x
x
x
=
+ + =
Krok 3.
3
2
1.4
x
x
h
= + =
2
3
3
3
3
( )
1
4.36
f x
x
x
=
+ + =
Krok 4.
4
3
1.6
x
x
h
= + =
2
4
4
4
4
(
)
1 5.16
f x
x
x
=
+ + =
Krok 5.
5
4
1.8
x
x
h
= + =
2
5
5
5
5
( )
1
6.04
f x
x
x
=
+ + =
Krok 6.
6
5
2
x
x
h
b
= + = =
2
6
6
6
6
( )
1
7
f x
x
x
=
+ + =
(
)
2
2
1
3
7
( )
0.2
3.64 4.36 5.16 6.04
4.84
2
2
2
2
n
a
b
i
i
f
f
f x dx
h
f
=
=
+
+
=
⋅
+
+
+
+
+
=
∑
∫
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
24
parabol
ę
– b
ę
dziemy oblicza
ć
sumy wycinków obszarów pod parabol
ą
. Algorytm jest
nast
ę
puj
ą
cy:
Rysunek 3. ZłoŜona metoda Simpsona.
Algorytm:
dzielimy przedział całkowania [a, b] na n równych odcinków;
n
a
b
h
−
=
,
(
)
(
) (
)
(
)
1
1
5
3
4
2
1
2
1
1
2
2
1
2
...
2
...
4
3
4
3
)
(
+
−
=
+
−
+
+
+
+
+
+
+
+
+
=
+
+
≅
∫
∑
n
n
n
b
a
n
i
i
i
i
f
f
f
f
f
f
f
f
h
f
f
f
h
dx
x
f
bł
ą
d dla metody Simpsona wyra
Ŝ
a si
ę
wzorem:
(
)
( )
( )
[ ]
b
a
f
h
a
b
R
,
,
180
1
4
4
∈
⋅
⋅
−
=
ς
ς
Przykład:
Oblicz warto
ść
całki
(
)
2
2
1
1
x
x
dx
+ +
∫
metod
ą
Simpsona, n = 6.
2
( )
1
f x
x
x
=
+ +
,
[ ]
1, 2
x
∈
,
n=6
Krok 1.
2 1
1
6
6
h
−
=
=
,
1
1
x
a
= =
2
1
1
1
1
( )
1
3
f x
x
x
=
+ + =
Krok 5.
5
4
5
3
x
x
h
= + =
2
5
5
5
5
49
( )
1
9
f x
x
x
=
+ + =
Krok 2.
2
1
7
6
x
x
h
= + =
2
2
2
2
2
127
(
)
1
36
f x
x
x
=
+ + =
Krok 6.
6
5
11
6
x
x
h
= + =
2
6
6
6
6
223
( )
1
36
f x
x
x
=
+ + =
25
Krok 3.
3
2
4
3
x
x
h
= + =
2
3
3
3
3
37
( )
1
9
f x
x
x
=
+ + =
Krok 7.
7
6
2
x
x
h
b
= + = =
2
7
7
7
7
(
)
1
7
f x
x
x
=
+ + =
Krok 4.
4
3
3
2
x
x
h
= + =
2
4
4
4
4
19
(
)
1
4
f x
x
x
=
+ + =
(
)
(
)
(
)
(
)
( )
( )
2
2
1
2
2
1
1
1
2
4
6
3
5
1
( )
4
3
4
2
3
1/ 6
127
19
223
37
49
29
3 4
2
7
4.8 3
3
36
4
36
9
9
6
n
b
i
i
i
i
a
n
h
f x dx
f
f
f
h
f
f
f
f
f
f
f
−
+
=
+
=
+ ⋅
+
=
+ ⋅
+ +
+ ⋅
+
+
=
+ ⋅
+
+
+ ⋅
+
+
=
=
∑
∫
12.
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,
cumtrapz
Całkowanie metod
ą
trapezów – metoda nieadaptacyjna
Q = quad(fun, a, b, [RelTol, AbsTol], trace)
fun
– funkcja podcałkowa,
a, b
– granice całkowania,
26
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)
13.
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
2
2
2
x
dx
+
∫
, n = 10
(
)
0
( )
Sin x
dx
π
∫
, n = 8
(
)
2
2
0
7
( ) 5
2
Cos x
x
dx
π
⋅
⋅
−
+
∫
, n = 10
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
Przykład implementacji funkcji w programie MATLAB – metoda prostokątów.
function
Int_MT = p_MT(fun, a, b, n)
if
abs(b - a) < eps | n <= 0,
Int_MT = 0;
return
;
end
h = (b - a)/n;
x = a + h/2:h:b-h/2;
fun_V = feval(fun, x);
Int_MT = h*sum(fun_V);
end
27
Przykład implementacji funkcji w programie MATLAB – metoda trapezów.
function
Int_MT = trapz_MT(fun, a, b, n)
if
abs(b - a) < eps | n <= 0,
Int_MT = 0;
return
;
end
h = (b - a)/n;
x = a + [0:n]*h;
fun_V = feval(fun, x);
Int_MT = h*((fun_V(1) + fun_V(n + 1))/2 + sum(fun_V(2:n)));
end
Przykład implementacji funkcji w programie MATLAB – metoda Simpsona.
function
Int_MS = simpson_MT(fun, a, b, n)
if
abs(b - a) < eps | n <= 0,
Int_MS = 0;
return
;
end
if
mod(n, 2) ~= 0,
n = n + 1;
end
h = (b - a)/n;
x = a + [0:n]*h;
fun_V = feval(fun, x);
fun_V(find(fun_V == inf)) = realmax;
fun_V(find(fun_V == -inf)) = -realmax;
w_odd = 2:2:n;
% w
ę
zły parzyste
w_even = 3:2:n - 1;
% w
ę
zły nieparzyste
Int_MS = h/3*(fun_V(1) + fun_V(n + 1) + 4*sum(fun_V(w_odd))
+ 2*sum(fun_V(w_even)));
end
[2]
Won Y. Yang, Wenwu Cao, Tae S. Chung, John Morris. Applied numerical methods
using MATLAB, John Wiley & Sons, Inc. 2005
28
KATEDRA URZĄDZEŃ I SYSTEMÓW AUTOMATYKI
LABORATORIUM
METOD NUMERYCZYCH
W TECHNICE
RÓWNANIA RÓśNICZKOWE
INSTRUKCJA LABORATORYJNA DO
Ć
WICZENIA NR 4.
WERSJA ROBOCZA
OPRACOWAŁ:
mgr in
Ŝ
. Tomasz KWA
Ś
NIEWSKI
POLITECHNIKA ŚWIĘTOKRZYSKA
KIELCE 2011
29
5.
Równanie róŜniczkowe
Ogólnie równanie ró
Ŝ
niczkowe zwyczajne pierwszego rz
ę
du mo
Ŝ
na zapisa
ć
:
( )
y
x
f
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
( )
( )
(
)
]
,
[
,
,
b
a
x
x
y
x
f
x
y
∈
∀
=
′
.
(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
( )
a
y
a
y
=
,
(3)
lub warunki brzegowe
( ) ( )
(
)
0
,
=
b
y
a
y
g
.
(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 znale
źć
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 znale
źć
rozwi
ą
zanie równania (1) w przedziale [x
0
, b]
spełniaj
ą
ce warunek pocz
ą
tkowy
(2)
y(x
0
) = y
0
Twierdzenie:
Je
ś
li funkcja f(x,y) jest ci
ą
gła w [x
0
, b] i spełnia w tym przedziale warunek
Lipschitza
(3)
( ) ( )
y
y
L
y
x
f
y
x
f
−
≤
−
~
,
~
,
30
Gdzie L- stała (stała Lipschitza), to w przedziale (x
0
, b) istaienje dokładnie jedno rozwi
ą
zanie
zagadnienia (1) i (2).
Całkuj
ą
c równanie (1) stronami, otrzymamy równanie całkowe
( )
( )
∫
+
=
b
x
dx
y
x
f
y
x
y
0
,
0
równowa
Ŝ
ne zagadnieniu (1) i (2).
Przykład metody Eulera
Rozwi
ą
za
ć
metod
ą
Eulera podany problem pocz
ą
tkowy:
2
'
2
y
x
y
= ⋅ +
,
( )
0
1
y
=
,
[
]
0, 0.5
x
∈
,
0.1
h
=
. Sprawd
ź
bł
ą
d metody dla rozwi
ą
zania
dokładnego
( )
2
4 5
4
2
x
y x
e
x
x
= − + ⋅ − ⋅ − ⋅
.
Wzór ogólny
metody Eulera
(
)
1
,
,
0,1, 2
i
i
i
i
y
y
h f x y
i
+
= + ⋅
=
K
KROK I.
( )
0
0
y x
y
=
, czyli
( )
0
1
y
=
,
0
0
x
=
,
0
i
=
(
)
(
)
(
)
2
2
1
0
0
0
0
0
0
,
2
1 0.1 2 0
1
1.1
y
y
h f x y
y
h
x
y
=
+ ⋅
=
+ ⋅ ⋅
+
= +
⋅ ⋅ + =
KROK II.
1
1.1
y
=
,
1
0
0.1
x
x
h
= + =
,
1
i
=
(
)
( )
(
)
2
2
1
1
1
,
1.1 0.1 2 0.1
1.1
1.212
y
y
h f x y
= + ⋅
=
+
⋅ ⋅
+
=
KROK III.
2
1.212
y
=
,
2
1
0.2
x
x
h
= + =
,
2
i
=
(
)
( )
(
)
2
3
2
2
2
,
1.212 0.1 2 0.2
1.212
1.3412
y
y
h f x y
=
+ ⋅
=
+
⋅ ⋅
+
=
KROK IV.
3
1.3412
y
=
,
3
2
0.3
x
x
h
= + =
,
3
i
=
(
)
( )
(
)
2
4
3
3
3
,
1.3412 0.1 2 0.3
1.3412
1.49332
y
y
h f x y
= + ⋅
=
+
⋅ ⋅
+
=
KROK V.
4
1.49332
y
=
,
4
3
0.4
x
x
h
= + =
,
4
i
=
(
)
( )
(
)
2
5
4
4
4
,
1.49332 0.1 2 0.4
1.49332
1.67465
y
y
h f x y
=
+ ⋅
=
+
⋅ ⋅
+
=
KROK VI.
5
1.67465
y
=
,
5
4
0.5
x
x
h
= + =
,
5
i
=
(
)
( )
(
)
2
6
5
5
5
,
1.67465 0.1 2 0.5
1.67465
1.89212
y
y
h f x y
= + ⋅
=
+
⋅ ⋅
+
=
31
Wyznaczamy bł
ą
d procentowy metody Eulera:
( )
2
4 5
4
2
x
y x
e
x
x
= − + ⋅ − ⋅ − ⋅
( )
0.5
0
0.660273
dok
y
y x dx
=
=
∫
( )
(
)
( )
(
)
(
)
(
)
6
0
0.660273 1
1.89212
100%
100% 13.96%
0.660273 1
0
dok
dok
y
y
y
y
y
δ
+
−
+ −
=
⋅
=
⋅
=
+
+
Rozwi
ą
za
ć
metod
ą
Runge-Kutty II rz
ę
du podany problem pocz
ą
tkowy:
2
'
2
y
x
y
= ⋅ +
,
( )
0
1
y
=
,
[
]
0, 0.5
x
∈
,
0.1
h
=
. Sprawd
ź
bł
ą
d metody dla rozwi
ą
zania
dokładnego
( )
2
4 5
4
2
x
y x
e
x
x
= − + ⋅ − ⋅ − ⋅
.
Wzór ogólny
metody Runge-Kutty II
rz
ę
du
KROK I.
( )
0
0
y x
y
=
, czyli
( )
0
1
y
=
,
0
0
x
=
,
0
i
=
(
)
(
)
(
)
(
)
(
) (
)
(
)
(
) (
)
(
)
2
2
1
0
0
0
0
2
2
2
0
0
1
0
0
1
,
2
0.1 2 0
1
0.1
,
2
0.1 2 0 0.1
1 0.1
0.112
g
h f x y
h
x
y
g
h f x
h y
g
h
x
h
y
g
= ⋅
= ⋅ ⋅
+
=
⋅ ⋅ + =
= ⋅
+
+
= ⋅ ⋅
+
+
+
=
⋅ ⋅ +
+ +
=
(
)
(
)
1
0
1
2
1
1
1
0.1 0.112
1.106
2
2
y
y
g
g
=
+ ⋅
+
= + ⋅
+
=
KROK II.
1
1.106
y
=
,
1
0
0.1
x
x
h
= + =
,
1
i
=
(
)
( )
(
)
(
)
(
) (
)
(
)
2
1
1
1
2
2
1
1
1
,
0.1 2 0.1
1.106
0.1126
,
0.1 2 0.1 0.1
1.106 0.1126
0.12986
g
h f x y
g
h f x
h y
g
= ⋅
=
⋅ ⋅
+
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
=
(
)
(
)
2
1
1
2
1
1
1.106
0.1126 0.12986
1.22723
2
2
y
y
g
g
= + ⋅
+
=
+ ⋅
+
=
KROK III.
2
1.22723
y
=
,
2
1
0.2
x
x
h
= + =
,
2
i
=
(
)
( )
(
)
(
)
(
) (
)
(
)
2
1
2
2
2
2
2
2
1
,
0.1 2 0.2
1.22723
0.130723
,
0.1 2 0.2 0.1
1.22723 0.130723
0.153795
g
h f x y
g
h f x
h y
g
= ⋅
=
⋅ ⋅
+
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
=
32
(
)
(
)
3
2
1
2
1
1
1.22723
0.130723 0.153795
1.36949
2
2
y
y
g
g
=
+ ⋅
+
=
+ ⋅
+
=
KROK IV.
3
1.36949
y
=
,
3
2
0.3
x
x
h
= + =
,
3
i
=
(
)
( )
(
)
(
)
(
) (
)
(
)
2
1
3
3
2
2
3
3
1
,
0.1 2 0.3
1.36949
0.154949
,
0.1 2 0.3 0.1
1.36949 0.154949
0.184
g
h f x y
g
h f x
h y
g
= ⋅
=
⋅ ⋅
+
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
=
(
)
(
)
4
3
1
2
1
1
1.36949
0.154949 0.184
1.53896
2
2
y
y
g
g
= + ⋅
+
=
+ ⋅
+
=
KROK V.
4
1.53896
y
=
,
4
3
0.4
x
x
h
= + =
,
4
i
=
(
)
( )
(
)
(
)
(
) (
)
(
)
2
1
4
4
2
2
4
4
1
,
0.1 2 0.4
1.53896
0.185896
,
0.1 2 0.4 0.1
1.53896 0.185896
0.222486
g
h f x y
g
h f x
h y
g
= ⋅
=
⋅ ⋅
+
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
=
(
)
(
)
5
4
1
2
1
1
1.53896
0.185896 0.222486
1.74315
2
2
y
y
g
g
=
+ ⋅
+
=
+ ⋅
+
=
Wyznaczamy bł
ą
d procentowy metody Runge-Kutty II rz
ę
du:
( )
2
4 5
4
2
x
y x
e
x
x
= − + ⋅ − ⋅ − ⋅
( )
0.5
0
0.660273
dok
y
y x dx
=
=
∫
( )
(
)
( )
(
)
(
)
(
)
5
0
0.660273 1
1.74315
100%
100%
4.99%
0.660273 1
0
dok
dok
y
y
y
y
y
δ
+
−
+ −
=
⋅
=
⋅
=
+
+
Rozwi
ą
za
ć
metod
ą
Runge-Kutty IV rz
ę
du podany problem pocz
ą
tkowy:
2
'
2
y
x
y
= ⋅ +
,
( )
0
1
y
=
,
[
]
0, 0.5
x
∈
,
0.1
h
=
. Sprawd
ź
bł
ą
d metody dla rozwi
ą
zania
dokładnego
( )
2
4 5
4
2
x
y x
e
x
x
= − + ⋅ − ⋅ − ⋅
.
Wzór ogólny
metody Runge-Kutty IV
rz
ę
du
KROK I.
( )
0
0
y x
y
=
, czyli
( )
0
1
y
=
,
0
0
x
=
,
0
i
=
33
(
)
(
)
(
)
2
2
1
0
0
0
0
2
2
1
1
2
0
0
0
0
2
2
2
3
0
0
0
0
,
2
0.1 2 0
1
0.1
0.1
0.1
,
2
0.1
2
0
1
0.1055
2
2
2
2
2
2
,
2
2
2
2
2
g
h f x y
h
x
y
g
g
h
h
g
h f x
y
h
x
y
g
g
h
h
g
h f x
y
h
x
y
= ⋅
= ⋅ ⋅
+
=
⋅ ⋅ + =
= ⋅
+
+
= ⋅ ⋅
+
+
+
=
⋅ ⋅ +
+ +
=
= ⋅
+
+
= ⋅ ⋅
+
+
+
(
)
(
) (
)
(
)
(
) (
)
(
)
2
2
2
4
0
0
3
0
0
3
0.1
0.1055
0.1 2
0
1
0.105775
2
2
,
2
0.1 2 0 0.1
1 0.105775
0.112578
g
h f x
h y
g
h
x
h
y
g
=
⋅ ⋅ +
+ +
=
= ⋅
+
+
= ⋅ ⋅
+
+
+
=
⋅ ⋅ +
+ +
=
(
)
(
)
1
0
1
2
3
4
1
1
2
2
1
0.1 2 0.1055 2 0.105775 0.112578
1.10585
6
6
y
y
g
g
g
g
=
+ ⋅
+ ⋅ + ⋅ +
= + ⋅
+ ⋅
+ ⋅
+
=
KROK II.
1
1.10585
y
=
,
1
0
0.1
x
x
h
= + =
,
1
i
=
(
)
( )
(
)
2
1
1
1
2
1
2
1
1
2
2
3
1
1
,
0.1 2 0.1
1.10585
0.112585
0.1
0.112585
,
0.1
2
0.1
1.10585
0.120714
2
2
2
2
0.1
0.120714
,
0.1
2
0.1
1.10585
2
2
2
2
g
h f x y
g
h
g
h f x
y
g
h
g
h f x
y
= ⋅
=
⋅ ⋅
+
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
(
)
(
) (
)
(
)
2
4
1
1
3
0.121121
,
0.1 2 0.1 0.1
1.10585 0.121121
0.130697
g
h f x
h y
g
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
=
(
)
(
)
2
1
1
2
3
4
1
1
2
2
1.10585
0.112585 2 0.120714 2 0.121121 0.130697
1.22701
6
6
y
y
g
g
g
g
= + ⋅
+ ⋅ + ⋅ +
=
+ ⋅
+ ⋅
+ ⋅
+
=
KROK III.
2
1.22701
y
=
,
2
1
0.2
x
x
h
= + =
,
2
i
=
(
)
( )
(
)
2
1
2
2
2
1
2
2
2
2
2
3
2
2
,
0.1 2 0.2
1.22701
0.130701
0.1
0.130701
,
0.1 2
0.2
1.22701
0.141736
2
2
2
2
0.1
0.141736
,
0.1
2
0.2
1.22701
2
2
2
2
g
h f x y
g
h
g
h f x
y
g
h
g
h f x
y
= ⋅
=
⋅ ⋅
+
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
(
)
(
) (
)
(
)
2
4
2
2
3
0.142288
,
0.1 2 0.2 0.1
1.22701 0.142288
0.15493
g
h f x
h y
g
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
=
(
)
(
)
3
2
1
2
3
4
1
1
2
2
1.22701
0.130701 2 0.141736 2 0.142288 0.15493
1.36929
6
6
y
y
g
g
g
g
= + ⋅
+ ⋅ + ⋅ +
=
+ ⋅
+ ⋅
+ ⋅
+
=
34
KROK IV.
3
1.36929
y
=
,
3
2
0.3
x
x
h
= + =
,
3
i
=
(
)
( )
(
)
2
1
3
3
2
1
2
3
3
2
2
3
3
3
,
0.1 2 0.3
1.36929
0.154929
0.1
0.154929
,
0.1
2
0.3
1.36929
0.169175
2
2
2
2
0.1
0.169175
,
0.1
2
0.3
1.36929
2
2
2
2
g
h f x y
g
h
g
h f
x
y
g
h
g
h f
x
y
= ⋅
=
⋅ ⋅
+
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
(
)
(
) (
)
(
)
2
4
3
3
3
0.169888
,
0.1 2 0.3 0.1
1.36929 0.169888
0.185918
g
h f x
h y
g
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
=
(
)
(
)
4
3
1
2
3
4
1
1
2
2
1.36929
0.154929 2 0.169175 2 0.169888 0.185918
1.53912
6
6
y
y
g
g
g
g
= + ⋅
+ ⋅ + ⋅ +
=
+ ⋅
+ ⋅
+ ⋅
+
=
KROK V.
4
1.53912
y
=
,
4
3
0.4
x
x
h
= + =
,
4
i
=
(
)
( )
(
)
2
1
4
4
2
1
2
4
4
2
2
3
4
4
,
0.1 2 0.4
1.53912
0.185912
0.1
0.185912
,
0.1
2
0.4
1.53912
0.203708
2
2
2
2
0.1
0.203708
,
0.1
2
0.4
1.53912
2
2
2
2
g
h f x y
g
h
g
h f
x
y
g
h
g
h f
x
y
= ⋅
=
⋅ ⋅
+
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
(
)
(
) (
)
(
)
2
4
4
4
3
0.204597
,
0.1 2 0.4 0.1
1.53912 0.204597
0.224372
g
h f x
h y
g
=
= ⋅
+
+
=
⋅ ⋅
+
+
+
=
(
)
(
)
5
4
1
2
3
4
1
1
2
2
1.53912
0.185912 2 0.203708 2 0.204597 0.224372
1.7436
6
6
y
y
g
g
g
g
=
+ ⋅
+ ⋅ + ⋅ +
=
+ ⋅
+ ⋅
+ ⋅
+
=
Wyznaczamy bł
ą
d procentowy metody Runge-Kutty IV rz
ę
du:
( )
2
4 5
4
2
x
y x
e
x
x
= − + ⋅ − ⋅ − ⋅
( )
0.5
0
0.660273
dok
y
y x dx
=
=
∫
35
( )
(
)
( )
(
)
(
)
(
)
5
0
0.660273 1
1.7436
100%
100%
5.01%
0.660273 1
0
dok
dok
y
y
y
y
y
δ
+
−
+ −
=
⋅
=
⋅
=
+
+
14.
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,
36
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:
2
'
2
y
x
y
= ⋅ +
,
( )
0
1
y
=
,
[
]
0, 0.5
x
∈
, 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:
(
)
t
y
y
y
⋅
=
⋅
−
⋅
+
10
sin
4
'
5
''
,
( )
( )
0
0
'
,
0
0
=
=
y
y
,
[ ]
3
,
0
∈
x
, 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
37
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:
2
'
2
y
x
y
= ⋅ +
,
( )
0
1
y
=
,
[
]
0, 0.5
x
∈
, h=0.01. Sprawd
ź
bł
ą
d metody dla rozwi
ą
zania
dokładnego
( )
2
4 5
4
2
x
y x
e
x
x
= − + ⋅ − ⋅ − ⋅
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
38
15.
Zadania
Zadanie 1.
Napisa
ć
funkcje w programie MATLAB rozwi
ą
zuj
ą
ce równania ró
Ŝ
niczkowe pierwszego
stopnia.
•
Metoda Runge-Kutty drugiego rz
ę
du
(
)
1
1
2
1
,
0,1, 2
2
i
i
y
y
g
g
i
+
= + ⋅
+
=
K
(
)
(
)
1
2
1
,
,
,
,
0,1, 2
i
i
i
i
g
h f x y
g
h f x
h y
g
i
= ⋅
= ⋅
+
+
=
K
•
Metoda Runge-Kutty czwartego rz
ę
du
(
)
1
1
2
3
4
1
2
2
,
0,1, 2
6
i
i
y
y
g
g
g
g
i
+
= + ⋅
+ ⋅ + ⋅ +
=
K
(
)
(
)
1
1
2
2
3
4
3
,
,
,
,
2
2
,
,
,
2
2
i
i
i
i
i
i
i
i
g
h
g
h f x y
g
h f
x
y
g
h
g
h f
x
y
g
h f x
h y
g
= ⋅
= ⋅
+
+
= ⋅
+
+
= ⋅
+
+
Zadanie 2.
Rozwi
ą
za
ć
równanie ró
Ŝ
niczkowe metodami poznanymi na zaj
ę
ciach.
•
'
5
3
y
x
y
= ⋅ + ⋅
,
( )
0
2
y
=
,
[ ]
0,1
x
∈
,
0.2
h
=
. Sprawd
ź
bł
ą
d metody dla rozwi
ą
zania
dokładnego
( )
(
)
3
1
5 23
15
9
x
y x
e
x
⋅
=
− + ⋅
− ⋅
.
•
3
'
y
x
y
= +
,
( )
0
1
y
=
,
[ ]
0,1
x
∈
,
0.2
h
=
. Sprawd
ź
bł
ą
d metody dla rozwi
ą
zania
dokładnego
( )
2
3
6 7
6
3
x
y x
e
x
x
x
= − + ⋅ − ⋅ − ⋅ −
.
•
2
'
2
y
x
y
=
+ ⋅
,
( )
0
2
y
=
,
[ ]
0,1
x
∈
,
0.2
h
=
. Sprawd
ź
bł
ą
d metody dla rozwi
ą
zania
dokładnego
( )
(
)
2
2
1
1 9
2
2
4
x
y x
e
x
x
⋅
=
− + ⋅
− ⋅ − ⋅
.
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.
39
KATEDRA URZĄDZEŃ I SYSTEMÓW AUTOMATYKI
LABORATORIUM
METOD NUMERYCZYCH
W TECHNICE
PIERWIASTKI RÓWNANIA
NIELINIOWEGO
INSTRUKCJA LABORATORYJNA DO
Ć
WICZENIA NR 5.
WERSJA ROBOCZA
OPRACOWAŁ:
mgr in
Ŝ
. Tomasz KWA
Ś
NIEWSKI
POLITECHNIKA ŚWIĘTOKRZYSKA
KIELCE 2011
40
Metoda równego podziału
(metoda bisekcji) – jedna z metod rozwi
ą
zywania równa
ń
nieliniowych. Opiera si
ę
ona na twierdzeniu Bolzano-Cauchy'ego:
Je
Ŝ
eli funkcja ci
ą
gła f(x) ma na ko
ń
cach przedziału domkni
ę
tego warto
ś
ci ró
Ŝ
nych znaków, to
wewn
ą
trz tego przedziału, istnieje co najmniej jeden pierwiastek równania f(x) = 0.
Aby mo
Ŝ
na było zastosowa
ć
metod
ę
równego podziału, musz
ą
by
ć
spełnione zało
Ŝ
enia:
•
Funkcja f(x) jest
ci
ą
gła
w przedziale domkni
ę
tym [a, b].
•
Funkcja przyjmuje ró
Ŝ
ne znaki na ko
ń
cach przedziału: f(a)*f(b) < 0.
Algorytm:
1.
Nale
Ŝ
y sprawdzi
ć
, czy pierwiastkiem równania jest punkt
2
i
a b
x
+
=
, czyli czy f(x
i
) = 0.
2.
Je
Ŝ
eli tak jest, algorytm ko
ń
czy si
ę
. W przeciwnym razie x
i
dzieli przedział [a, b] na dwa
mniejsze przedziały [a, x
i
] i [x
i
, b].
3.
Nast
ę
pnie wybierany jest ten przedział, dla którego spełnione jest drugie zało
Ŝ
enie, tzn.
albo f(a)*f(x
i
) < 0 albo f(x
i
)*f(b) < 0. Cały proces powtarzany jest dla wybranego przedziału.
Algorytm ko
ń
czy si
ę
w punkcie 2, lub jest przerywany, gdy zaw
ęŜ
ony przedział b
ę
dzie
mniejszy od zało
Ŝ
onej tolerancji, tzn. gdy
i
i
b
a
−
≤ ∆
.
Przykład:
Znale
źć
rozwi
ą
zanie równania
2
( )
2
2
2
f x
x
x
= ⋅ + ⋅ −
w przedziale
[ ] [
]
,
-2, 2
a b
=
metod
ą
bisekcji
zakładaj
ą
c wielko
ść
bł
ę
du
0.25
∆ =
.
KROK I.
Sprawdzamy:
b a
− ≤ ∆
, b=2, a=-2, otrzymujemy
( )
2
2
4
0.25
− − = > ∆ =
.
Wyznaczmy
1
2
2
0
2
2
a b
x
+
− +
=
=
=
, nast
ę
pnie sprawdzamy czy
( )
1
0
f x
=
, otrzymujemy
( )
( )
2
0
2 0
2 0 2
2
f
= ⋅
+ ⋅ − = −
czyli
( )
1
0
f x
≠
.
Kolejny krok to okre
ś
lenie nowego przedziału
[ ] [
]
,
-2, 2
a b
=
, takiego dla którego spełniony
jest warunek
( ) ( )
0
f a
f b
⋅
<
. Powstaj
ą
dwa nowe przedziały w postaci:
[ ] [
]
1
,
-2, 0
a x
=
, sprawdzamy
( ) ( )
( )
2
0
2
2
4
0
f
f
− ⋅
= ⋅ − = − <
,
[ ] [ ]
1
,
0, 2
x b
=
, sprawdzamy
( ) ( )
0
2
2 10
20
0
f
f
⋅
= − ⋅ = − <
.
Warunek
( ) ( )
0
f a
f b
⋅
<
spełniony jest dla obu przedziałów oznacza to,
Ŝ
e funkcja f(x)
posiada miejsca zerowe w obu przedziałach. W pierwszej kolejno
ś
ci sprawdzimy przedział
[ ] [
]
1
,
-2, 0
a x
=
tym samym otrzymujemy nowy przedział
[ ] [
]
,
-2, 0
a b
=
i powtarzamy cały
algorytm oblicze
ń
dla nowego przedziału.
41
KROK II.
Sprawdzamy:
b a
− ≤ ∆
, b=0, a=-2, otrzymujemy
( )
0
2
2
0.25
− − = > ∆ =
.
Wyznaczmy
1
2 0
1
2
2
a b
x
+
− +
=
=
= −
, nast
ę
pnie sprawdzamy czy
( )
1
0
f x
=
, otrzymujemy
( )
( )
( )
2
1
2
1
2
1
2
2
f
− = ⋅ −
+ ⋅ − − = −
czyli
( )
1
0
f x
≠
.
Kolejny krok to okre
ś
lenie nowego przedziału
[ ] [
]
,
-2, 0
a b
=
, takiego dla którego spełniony
jest warunek
( ) ( )
0
f a
f b
⋅
<
. Powstaj
ą
dwa nowe przedziały w postaci:
[ ] [
]
1
,
-2, -1
a x
=
, sprawdzamy
( ) ( )
( )
2
1
2
2
4
0
f
f
− ⋅
− = ⋅ − = − <
,
[ ] [
]
1
,
-1, 0
x b
=
, sprawdzamy
( ) ( )
( )
1
0
2
2
4
0
f
f
− ⋅
= − ⋅ − = >
.
Warunek
( ) ( )
0
f a
f b
⋅
<
spełniony jest dla przedziału
[ ] [
]
1
,
-2, -1
a x
=
oznacza to,
Ŝ
e funkcja
f(x) posiada miejsca zerowe w tym przedziale. Wybieramy przedział
[ ] [
]
1
,
-2, -1
a x
=
tym
samym otrzymujemy nowy przedział
[ ] [
]
,
-2, -1
a b
=
i powtarzamy cały algorytm oblicze
ń
dla
nowego przedziału.
KROK III.
Sprawdzamy:
b a
− ≤ ∆
, b=-1, a=-2, otrzymujemy
( )
1
2
1
0.25
− − − = > ∆ =
.
Wyznaczmy
( )
1
2
1
1.5
2
2
a b
x
− + −
+
=
=
= −
,
nast
ę
pnie
sprawdzamy
czy
( )
1
0
f x
=
,
otrzymujemy
(
)
(
)
(
)
2
1.5
2
1.5
2
1.5
2
0.5
f
−
= ⋅ −
+ ⋅ −
− = −
czyli
( )
1
0
f x
≠
.
Kolejny krok to okre
ś
lenie nowego przedziału
[ ] [
]
,
-2, -1.5
a b
=
, takiego dla którego
spełniony jest warunek
( ) ( )
0
f a
f b
⋅
<
. Powstaj
ą
dwa nowe przedziały w postaci:
[ ] [
]
1
,
-2, -1.5
a x
=
, sprawdzamy
( ) (
)
(
)
2
1.5
2
0.5
1 0
f
f
− ⋅
−
= ⋅ −
= − <
,
[ ] [
]
1
,
-1.5, -1
x b
=
, sprawdzamy
(
) ( )
( )
1.5
1
0.5
2
1 0
f
f
−
⋅
− = −
⋅ − = >
.
Warunek
( ) ( )
0
f a
f b
⋅
<
spełniony jest dla przedziału
[ ] [
]
1
,
-2, -1.5
a x
=
oznacza to,
Ŝ
e
funkcja f(x) posiada miejsca zerowe w tym przedziale. Wybieramy przedział
[ ] [
]
1
,
-2, -1.5
a x
=
tym samym otrzymujemy nowy przedział
[ ] [
]
,
-2, -1.5
a b
=
i powtarzamy
cały algorytm oblicze
ń
dla nowego przedziału.
KROK IV.
Sprawdzamy:
b a
− ≤ ∆
, b=-1.5, a=-2, otrzymujemy
( )
1.5
2
0.5
0.25
−
− − =
> ∆ =
.
Wyznaczmy
(
)
1
2
1.5
1.75
2
2
a b
x
− + −
+
=
=
= −
, nast
ę
pnie sprawdzamy czy
( )
1
0
f x
=
,
otrzymujemy
(
)
(
)
(
)
2
1.75
2
1.75
2
1.75
2
0.625
f
−
= ⋅ −
+ ⋅ −
− =
czyli
( )
1
0
f x
≠
.
Kolejny krok to okre
ś
lenie nowego przedziału
[ ] [
]
,
-2, -1.75
a b
=
, takiego dla którego
spełniony jest warunek
( ) ( )
0
f a
f b
⋅
<
. Powstaj
ą
dwa nowe przedziały w postaci:
[ ] [
]
1
,
-2, -1.75
a x
=
, sprawdzamy
( ) (
)
2
1.75
2 0.625 1.25
0
f
f
− ⋅
−
= ⋅
=
>
,
[ ] [
]
1
,
-1.75, -1.5
x b
=
, sprawdzamy
(
) (
)
(
)
1.75
1.5
0.625
0.5
0.3125
0
f
f
−
⋅
−
=
⋅ −
= −
<
.
42
Warunek
( ) ( )
0
f a
f b
⋅
<
spełniony jest dla przedziału
[ ] [
]
1
,
-1.75, -1.5
x b
=
oznacza to,
Ŝ
e
funkcja f(x) posiada miejsca zerowe w tym przedziale. Wybieramy przedział
[ ] [
]
1
,
-1.75, -1.5
x b
=
tym samym otrzymujemy nowy przedział
[ ] [
]
,
-1.75, -1.5
a b
=
i
powtarzamy cały algorytm oblicze
ń
dla nowego przedziału.
KROK V.
Sprawdzamy:
b a
− ≤ ∆
, b=-1.5, a=-1.75, otrzymujemy
(
)
1.5
1.75
0.25
0.25
−
− −
=
= ∆ =
.
Wyznaczmy
(
)
1
1.75
1.5
1.625
2
2
a b
x
−
+ −
+
=
=
= −
. Miejscem zerowym funkcji f(x) jest punkt
*
1.625
x
= −
wyznaczony z dokładno
ś
ci
ą
0.25
∆ =
.
Kolejny krok to powtórzenie całego algorytmu dla przedziału
[ ] [ ]
1
,
0, 2
x b
=
z kroku
pierwszego.
PRZYKŁAD
Przy pomocy metody bisekcji wyznaczy
ć
pierwiastki równania:
( )
1
6
5
2
+
⋅
−
⋅
−
=
x
x
x
f
function
[x, err, xx] = bisct(f, a, b, TolX, MaxIter)
fa = feval(f, a);
fb = feval(f, b);
for
k = 1: MaxIter
xx(k) = (a + b)/2;
fx = feval(f, xx(k));
err = (b-a)/2;
if
abs(fx) < eps | abs(err)<TolX,
break
;
elseif
fx*fa > 0,
a = xx(k);
fa = fx;
else
b = xx(k);
end
end
x = xx(k);
if
k == MaxIter,
fprintf(
'Najlepszy wynik w %d iteracji\n'
, MaxIter),
end
- skrypt z wywołaniem funkcji bisekcji
function
test_bisct
a = 0;
% dolna granica
b = 2;
% górna granica
tol = 1e-1;
%
Ŝą
dana dokładno
ść
iter = 10;
% ilo
ść
iteracji funkcji 'bisct'
fun = @(x)(-5*x.^2 -6*x +1);
% równanie funkcji
[x, delta, x_val] = bisct(fun, a, b, tol, iter)
% wywołanie funkcji
end
43
Metoda Newtona
(metoda stycznych) –
algorytm
wyznaczania przybli
Ŝ
onej warto
ś
ci
pierwiastka funkcji
y = f(x) jednej zmiennej w zadanym
przedziale
[a, b]. Zało
Ŝ
enia metody
s
ą
nast
ę
puj
ą
ce:
•
W przedziale [a, b] znajduje si
ę
dokładnie jeden pierwiastek.
•
Funkcja ma ró
Ŝ
ne znaki na kra
ń
cach przedziału, f(a)*f(b) < 0.
•
Pierwsza i druga
pochodna
maj
ą
stały znak w tym przedziale.
W pierwszym kroku metody wybierany jest ten kraniec przedziału, dla którego znak funkcji i
drugiej pochodnej s
ą
równe, a nast
ę
pnie z tego punktu (albo (a, f(a)) albo (b, f(b)))
wyprowadzana jest
styczna
.
Odci
ę
ta
punktu przeci
ę
cia stycznej z osi
ą
OX jest pierwszym
przybli
Ŝ
eniem rozwi
ą
zania. Je
ś
li to przybli
Ŝ
enie nie jest satysfakcjonuj
ą
ce, wówczas punkt
(x
1
, f(x
1
)) jest wybierany jako koniec przedziału i wszystkie czynno
ś
ci s
ą
powtarzane. Proces
jest kontynuowany, a
Ŝ
zostanie uzyskane wystarczaj
ą
co dobre przybli
Ŝ
enie pierwiastka
(warto
ść
funkcji w wyznaczonym punkcie).
Przykład:
Znale
źć
rozwi
ą
zanie równania
2
( )
2
2
2
f x
x
x
= ⋅ + ⋅ −
w przedziale
[ ] [
]
,
-2, 2
a b
=
metod
ą
stycznych
zakładaj
ą
c wielko
ść
bł
ę
du
0.25
∆ =
.
KROK I.
Wybieramy punkt startowy z warunku
( ) ( )
"
0
0
0
f x
f
x
⋅
≥
podstawiamy kolejno:
( ) ( )
"
0
2,
2
2
2 4
8
0
x
a
f
f
= = −
− ⋅
− = ⋅ = >
,
( ) ( )
"
0
2,
2
2
10 4
40
0
x
b
f
f
= =
⋅
= ⋅ =
>
.
Oba warunki s
ą
spełnione wobec tego z obu punktów startowych otrzymamy miejsce zerowe
funkcji f(x) w tym przypadku wybieramy dowolny z punktów startowych, przyjmijmy,
Ŝ
e
punkt startowy
0
2
x
b
= =
.
Po pierwsze obliczamy warto
ś
ci
( )
0
f x
i
( )
'
0
f
x
:
( )
( )
2
0
2 2
2 2 2 10
f x
= ⋅
+ ⋅ − =
,
( )
'
0
4 2 2
8
f
x
= ⋅ + =
.
Nast
ę
pnie obliczamy nowe przybli
Ŝ
enie rozwi
ą
zania wykorzystuj
ą
c wzór metody stycznych:
( )
( )
0
1
0
'
0
10
2
0.75
8
f x
x
x
f
x
= −
= −
=
.
Kolejny podpunkt to sprawdzenie czy
0
1
x
x
−
≤ ∆
. Czyli
2 0.75
1.25
0.25
−
=
> ∆ =
.
Warunek nie jest spełniony, powstaje nowy przedział
[ ] [
]
,
-2, 0.75
a b
=
.
KROK II.
Wybieramy punkt startowy z warunku
( ) ( )
"
0
0
0
f x
f
x
⋅
≥
podstawiamy kolejno:
( ) ( )
"
0
2,
2
2
2 4
8
0
x
a
f
f
= = −
− ⋅
− = ⋅ = >
,
(
) (
)
"
0
0.75,
0.75
0.75
0.625 5
3.125
0
x
b
f
f
= =
⋅
=
⋅ =
>
.
Oba warunki s
ą
spełnione wobec tego z obu punktów startowych otrzymamy miejsce zerowe
funkcji f(x) w tym przypadku wybieramy dowolny z punktów startowych, przyjmijmy,
Ŝ
e
punkt startowy
0
0.75
x
b
= =
.
Po pierwsze obliczamy warto
ś
ci
( )
0
f x
i
( )
'
0
f
x
:
( )
(
)
2
0
2 0.75
2 0.75 2
0.625
f x
= ⋅
+ ⋅
− =
,
( )
'
0
4 0.75 2
5
f
x
= ⋅
+ =
.
Nast
ę
pnie obliczamy nowe przybli
Ŝ
enie rozwi
ą
zania wykorzystuj
ą
c wzór metody stycznych:
44
( )
( )
0
1
0
'
0
0.625
0.75
0.625
5
f x
x
x
f
x
= −
=
−
=
.
Kolejny podpunkt to sprawdzenie czy
0
1
x
x
−
≤ ∆
. Czyli
0.75 0.625
0.125
0.25
−
=
< ∆ =
.
Wyznaczony punkt
*
0.625
x
=
jest miejscem zerowym funkcji f(x) z dokładno
ś
ci
ą
0.25
∆ =
.
Kolejny krok to powtórzenie całego algorytmu dla punktu startowego
0
2
x
a
= = −
.
PRZYKŁAD
- skrypt z wywołaniem funkcji newtona
function
test_newton
x0 = -5;
% punkt startowy
tol = 1e-3;
%
Ŝą
dana dokładno
ść
iter = 10;
% ilo
ść
iteracji funkcji 'newton'
fun = @(x)(- 5*x.^2 - 6*x + 1);
% równanie funkcji
dfun = @(x)(- 10*x - 6);
% równanie pochodnej funkcji
[x, delta, x_val] = newton(fun, dfun, x0, tol, iter)
% wywołanie
funkcji
end
Zadania
Zadanie 1.
Napisa
ć
funkcje w programie MATLAB rozwi
ą
zuj
ą
ce metody poszukiwania pierwiastków
równa
ń
nieliniowych przedstawione w instrukcji.
•
sprawdzi
ć
działanie funkcji ‘
bisct
’,
•
poprawi
ć
funkcj
ę
‘
bisct
’ według wskazówek prowadz
ą
cego,
•
napisa
ć
funkcj
ę
w programie MATLAB metody stycznych.
Zadanie 2.
Wyznaczy
ć
pierwiastki równania metodami poznanymi w instrukcji.
•
2
( )
2
4
9
f x
x
x
= − ⋅ + ⋅ +
,
[
]
-2, 6
x
∈
,
0.2
∆ =
.
•
4
3
2
1
( )
2
5
4
f x
x
x
x
= −
− ⋅ + −
,
[
]
-10, 4
x
∈
,
0.4
∆ =
.
•
2
x
4
)
x
(
Sin
)
x
(
f
−
+
=
,
[
]
-2, 2
x
∈
,
0.1
∆ =
.
Narysuj wykres funkcji i zaznacz na wykresie miejsce zerowe funkcji.
Zadanie 3.
Korzystaj
ą
z metod wbudowanej w programie MATLAB ‘fsolve’ nale
Ŝ
y przygotowa
ć
skrypt
do wyznaczania pierwiastków równa
ń
z zadania 2.
45
KATEDRA URZĄDZEŃ I SYSTEMÓW AUTOMATYKI
LABORATORIUM
METOD NUMERYCZYCH
W TECHNICE
UKŁAD RÓWNAŃ NIELINIOWYCH
INSTRUKCJA LABORATORYJNA DO
Ć
WICZENIA NR 6.
WERSJA ROBOCZA
OPRACOWAŁ:
mgr in
Ŝ
. Tomasz KWA
Ś
NIEWSKI
POLITECHNIKA ŚWIĘTOKRZYSKA
KIELCE 2011
46
Układ równań nieliniowych
Uogólnienie metody Newtona-Raphsona na przypadek wielowymiarowy to metoda Newtona,
któr
ą
mo
Ŝ
na zapisa
ć
jako:
( ) ( )
1
1
,
0,1, 2
i
i
i
i
x
x
x
f x
i
+
−
= −
⋅
=
J
K
gdzie:
( )
( ) ( )
( )
1
2
,
,
,
i
i
i
i
T
n
f x
f x
f
x
f
x
=
K
- układ równa
ń
zapisany w postaci wektora
( )
1
1
1
1
1
i
n
n
n
n
f
f
x
x
x
f
f
x
x
−
∂
∂
∂
∂
=
∂
∂
∂
∂
J
K
K
K
K
K
- macierz Jacobiego,
( )
(
)
1
det
0
i
x
−
≠
J
[
]
*
1
2
,
,
,
T
n
x
x x
x
=
K
- szukane zmienne niezale
Ŝ
ne
Przykład metody Newtona
Rozwi
ą
za
ć
metod
ą
Newtona podany układ równa
ń
dwóch zmiennych:
(
)
(
)
2
3
1
1
2
1
1
2
2
3
2
2
1
2
1
1
2
2
,
4
5
1.5
,
3
3.5
f x x
x
x
x
x
f
x x
x
x x
x
=
− ⋅ ⋅
+ ⋅
=
=
+ ⋅ ⋅ −
=
(
)
⋅
−
⋅
⋅
+
⋅
⋅
+
⋅
−
−
⋅
=
2
1
2
2
1
2
2
2
1
2
1
2
1
2
3
3
3
15
2
4
4
2
,
x
x
x
x
x
x
x
x
x
x
x
J
KROK I.
0
0
1
2
1,
1
x
x
=
=
- zerowe przybli
Ŝ
enie,
0
i
=
(
)
(
)
0
0
1
0
0
1
2
1
2
2 13
1
13
1
,
,
,
6
1
6
2
80
x x
x x
−
−
−
=
= −
⋅
−
−
J
J
47
(
) (
)
(
)
0
0
1
0
1
1
2
1
1
1
0
0
1
2
1
0
0
0
2
2
2
1
2
,
,
,
1
1
13
0.5
1.0875
1
1
6
2
0.5
0.975
80
f x x
x
x
x x
x
x
f
x x
−
=
−
⋅
=
−
− −
⋅
⋅
=
−
−
−
J
KROK II.
1
1
1
2
1.0875,
0.975
x
x
=
=
,
1
i
=
(
) (
)
(
)
1
1
2
1
1
1
2
1
1
1
1
1
1
2
2
1
1
1
2
2
2
1
2
,
,
,
208
662
1.0875
0.0216
1.085386
12737
4413
0.975
545
66
0.0164
0.972891
6767
2989
f
x x
x
x
x x
x
x
f
x x
−
=
−
⋅
=
−
−
⋅
=
J
KROK III.
2
2
1
2
1.085386,
0.972891
x
x
=
=
,
2
i
=
(
) (
)
(
)
2
2
3
2
1
1
2
1
1
1
2
2
1
2
3
2
2
2
2
2
2
1
2
,
,
,
231
651
1.085386
0.0000594
1.085383
14057
4327
0.972891
403
389
0.0000227
0.972885
4980
17479
f x x
x
x
x x
x
x
f
x x
−
=
−
⋅
=
−
−
⋅
=
J
Odpowied
ź
:
[
]
*
1.085383, 0.972885
T
x
=
(
)
(
)
3
3
2
3
1
1
2
1
1
2
2
3
3
3
2
2
1
2
1
1
2
2
,
4
5
1.5
4.0765
10
,
3
3.5
3.0261
11
f x x
x
x
x
x
e
f
x x
x
x x
x
e
=
− ⋅ ⋅
+ ⋅
−
=
⋅ −
=
+ ⋅ ⋅ −
−
=
⋅ −
48
Implementacja metody Newtona w programie MATLAB.
•
funkcja Newtona
function
[x, fx, xx] = newton_urn(fun, x0, TolX, MaxIter)
%input:
% f - układ równa
ń
zapisany w postaci wektora
% x0 - wektor punktów startowych dla zmiennych
% TolX - the upper limit of |x(k) - x(k - 1)|
% MaxIter - maksymalna liczba iteracji
%output:
% x - warto
ś
ci zmiennych obliczonych
% fx - warto
ś
ci równa
ń
dla zmiennych obliczonych
% xx - "historia" zmiennych obliczonych
h = 1e-4; TolFun = eps; EPS = 1e-6;
fx = feval(fun, x0);
Liczba_f = length(fx); Liczba_x = length(x0);
% sprawdzenie czy ilo
ść
równa
ń
jest równa ilo
ś
ci punktów startowych
if
Liczba_f ~= Liczba_x,
error(
'Niepoprawne dane! Liczba równa
ń
nie jest równa ilo
ś
ci punktów
startowych'
);
end
% je
Ŝ
eli ilo
ść
parametrów wej
ś
ciowych funkcji jest mniejsza od 4 to:
if
nargin < 4,
MaxIter = 100;
end
% je
Ŝ
eli ilo
ść
parametrów wej
ś
ciowych funkcji jest mniejsza od 3 to:
if
nargin < 3,
TolX = EPS;
end
xx(1, :) = x0(:).';
% przygotowanie wektora k-tego kroku
for
k = 1: MaxIter
dx = -jacobian(fun, xx(k,:), h)\fx(:);
xx(k + 1,:) = xx(k,:) + dx.';
fx = feval(fun, xx(k + 1,:));
fxn = norm(fx);
if
fxn < TolFun | norm(dx) < TolX,
break
;
end
end
x = xx(k + 1,:);
% wektor szukanych zminnych
if
k == MaxIter,
fprintf(
'Wynik otrzymany w %d iteracji\n'
, MaxIter),
end
49
•
funkcja obliczaj
ą
ca Jacobian
function
g = jacobian(fun, x, h)
if
nargin < 3,
h = 1e-3;
end
h2 = 2*h;
N = length(x);
x = x(:).';
I = eye(N);
for
n = 1:N
g(:, n) = (feval(fun, x + I(n, :)*h) - feval(fun, x - I(n, :)*h))'/h2;
end
•
skrypt testuj
ą
cy funkcj
ę
Newtona
Rozwi
ą
za
ć
przy pomocy funkcji
„
newton_urn
” podany układ równa
ń
dwóch zmiennych:
(
)
(
)
2
3
1
1
2
1
1
2
2
3
2
2
1
2
1
1
2
2
,
4
5
1.5
,
3
3.5
f x x
x
x
x
x
f
x x
x
x x
x
=
− ⋅ ⋅
+ ⋅
=
=
+ ⋅ ⋅ −
=
function
test_newton_urn
x0 = [2 3];
[x, f_x, x_x] = newton_urn(@uklad_rownan, x0)
function
y = uklad_rownan(x)
y(1) = x(1).^2 - 4*x(1).*sqrt(x(2)) + 5*x(2).^3 - 1.5;
y(2) = x(1).^3 + 3*x(1).*x(2) - x(2).^2 - 3.5;
end
end
Zadania
Zadanie 1.
Rozwi
ąŜ
układ równa
ń
dwóch zmiennych metod
ą
Newtona dla nast
ę
puj
ą
cych punktów
startowych:
a -------------
0
0
1
2
2,
1
x
x
=
= −
(
)
(
)
5
3
1
1
2
1
1
1
2
3
2
1
2
2
2
1
,
5
8.4
8
10
,
16
8
4
f x x
x
x
x
x
f
x x
x
x
x
= ⋅
−
⋅
+ ⋅ + =
= ⋅
− ⋅ + =
50
b -------------
0
0
1
2
4,
5
x
x
=
=
(
)
(
)
2
2
1
1
2
1
2
1
2
2
2
2
1
2
1
2
1
,
2
48
40
304
,
3
69
f x x
x
x
x
x
f
x x
x
x
x
= ⋅
+
− ⋅ − ⋅ = −
= − − ⋅
+ = −
c -------------
0
0
1
2
2,
1
x
x
=
= −
(
)
(
)
3
2
1
1
2
1
2
1
2
2
2
2
1
2
2
1
2
1
,
4
2
7
3
,
3
6
6
f
x x
x
x
x x
x
f
x x
x
x x
x
= − ⋅
+ ⋅
+ ⋅ − ⋅
=
= ⋅
+ ⋅ ⋅ − =
Zadanie 2.
Korzystaj
ą
z metod wbudowanej w programie MATLAB funkcji ‘fsolve’ nale
Ŝ
y przygotowa
ć
skrypt do wyznaczania pierwiastków równa
ń
z zadania 1.