Analiza matematyczna - elementy obliczeń

(1)

Wydział Matematyki i Nauk Informacyjnych Politechnika Warszawska

Analiza matematyczna - elementy obliczeń

Janusz Wąsowski

(2)

Analiza matematyczna, II semestr.

Elementy oblicze´n

Janusz W ˛asowski

1 Przybli˙zone obliczanie całek

Niech f : [a, b] → R b ˛edzie funkcj ˛a całkowaln ˛a (w sensie Riemanna) i niech I(f) oznacza całk ˛e oznaczon ˛a Riemanna funkcji f na [a, b], tzn.

I(f) =

b a

f(x)dx.

W niniejszym punkcie rozpatrzymy zagadnienie przybli˙zonego obliczania takich całek.

W kursie analizy matematycznej omawiane s ˛a ró˙zne analityczne sposoby obliczania takich całek.

Polegaj ˛a one na znajdowaniu funkcji pierwotnej funkcji podcałkowej poprzez wykorzystanie ró˙znych przekształceń całki, głównie przekształceń opartych na wzorach całkowania przez cz ˛e´sci i przez podstawienie. Sposoby te daj ˛a si ˛e zastosować do do´sć w ˛askiej klasy całek; poza t ˛a klas ˛a zachodzi konieczno´sć obliczania przybli˙zonej warto´sci całki metodami numerycznymi. Ponadto, je˙zeli funkcja podcałkowa jest okre´slona za pomoc ˛a tablicy, to mo˙zemy obliczyć tylko przybli˙zon ˛a warto´sć całki.

1.1 Sumy całkowe: interpretacja geometryczna, własno´sci Oznaczenia stosowane w definicji całki.

• ∆ = (x0, x1, ..., xn) - podział przedziału [a, b], przy czym a= x0 < x1 < ... < xn−1< xn= b;

• d(∆) = max {|xi− xi−1| | i = 1, 2, ..., n} - ´srednica podziału ∆;

• X∆= (ξ₁, ξ₂, ..., ξ_n) , gdzie ξ_ijest punktem po´srednim, wybranym z ka˙zdego z przedziałów [xi−1, xi] podziału ∆;

(3)

• suma całkowa funkcji f dla podziału ∆ i odpowiadaj ˛acemu mu wyborowi punktów po´srednich X∆:

σ(f, ∆, X∆) =

n i=1

f(ξ_i)(xi− xi−1);

• suma dolna i suma górna funkcji f dla podziału ∆ : s(f, ∆) =

n i=1

mi(xi− xi−1), S(f, ∆) =

n i=1

Mi(xi− xi−1),

gdzie mi= inf

x∈[x_i−1,xi]f(x), Mi = sup

x∈[x_i−1,xi]

f(x).

Interpretacja geometryczna.

Suma całkowa σ(f, ∆, X∆), gdy f ≥ 0, jest w interpretacji geometrycznej sum ˛a pól prostok ˛atów o podstawach (xi− xi−1) i wysoko´sciach f(ξ_i). Na rysunku podano interpretacj ˛e geometryczn ˛a sumy całkowej w przypadku n = 3.

Pole P obszaru domkni ˛etego ograniczonego wykresem funkcji y = f (x) ≥ 0, osi ˛a Ox i prostymi x= a, x = b jest granic ˛a sumy pól takich prostok ˛atów, gdy ´srednica podziału d(∆) → 0

P = lim

d(∆)→0

n i=1

f(ξ_i)(xi− xi−1) = I(f).

Ró˙znica sum S(f, ∆)−s(f, ∆) jest sum ˛a pól prostok ˛atów o podstawach (xi−xi−1) i wysoko´sciach Mi− mi :

S(f, ∆) − s(f, ∆) =

n i=1

(Mi− mi)(xi− xi−1). (1) Na rysunku podano interpretacj ˛e geometryczn ˛a tej ró˙znicy w przypadku n = 3.

(4)

Sumy całkowe - ocena bł ˛edu.

Mamy zapewnion ˛a zbie˙zno´s´c

σ(f, ∆, X∆) →

d(∆)→0I(f ), (2)

wi ˛ec σ(f, ∆, X∆) przybli˙za całk ˛e I(f) z dowoln ˛a dokładno´sci ˛a, je´sli tylko ´srednica podziału d(∆) jest wystarczaj ˛aco mała. Obecnie postawimy pytanie o szybko´s´c tej zbie˙zno´sci. Zale˙zy ona zarówno od podziału ∆, wyboru punktów po´srednich X∆, jak i regularno´sci funkcji f.

Ocenimy zbie˙zno´s´c w (2) szacuj ˛ac bł ˛ad

|I(f ) − σ(f, ∆, X∆)| (3)

w zale˙zno´sci od d(∆). Nasze oceny oprzemy na nast ˛epuj ˛acym znanym oszacowaniu

|I(f ) − σ(f, ∆, X∆)| ≤ S(f, ∆) − s(f, ∆). (4) Zauwa˙zmy, ˙ze ró˙znica S(f, ∆)−s(f, ∆) nie zale˙zy od X∆(nie zale˙zy od sposobu wyboru punktów po´srednich ξ_i). Ponadto wiemy o zbie˙zno´sci

S(f, ∆) − s(f, ∆) −→

d(∆)→00. (5)

Obecnie ocenimy t ˛e zbie˙zno´s´c w zale˙zno´sci od d(∆), gdy f jest funkcj ˛a monotoniczn ˛a lub funkcj ˛a ci ˛agł ˛a. Niech f b ˛edzie funkcj ˛a monotoniczn ˛a. Korzystaj ˛ac z tego, ˙ze

Mi− mi = |f(xi) − f (xi−1)|

mamy

S(f, ∆) − s(f, ∆) =

n i=1

|f(xi) − f (xi−1)| · (xi− xi−1) ≤

≤ d(∆)

n i=1

|f(xi) − f (xi−1)| = d(∆) · |f(b) − f(a)| ,

(5)

a st ˛ad otrzymujemy nast ˛epuj ˛ace oszacowanie bł ˛edu bezwzgl ˛ednego przybli˙zenia σ(f, ∆, X∆)

|I(f) − σ(f, ∆, X∆)| ≤ |f(b) − f (a)| · d(∆). (6) Do oceny zbie˙zno´sci w (5), gdy f jest funkcj ˛a ci ˛agł ˛a, u˙zyteczny jest moduł ci ˛agło´sci funkcji ωf. Wykonamy szacowania

Mi− mi = sup

x∈[xi−1,xi]

f(x) − inf

x∈[xi−1,xi]f(x) =

= sup

x,y∈[x_i−1,xi]

|f(x) − f(y)| ≤ ωf(xi− xi−1) ≤ ωf(d(∆)).

Opieraj ˛ac si ˛e na przedstawieniu (1), mamy S(f, ∆) − s(f, ∆) ≤ ωf(d(∆))

n i=1

(xi− xi−1) ≤ (b − a) · ωf(d(∆)). (7)

Ostatecznie, z (4) i (7) wynika nast ˛epuj ˛ace oszacowanie bł ˛edu

|I(f) − σ(f, ∆, X∆)| ≤ (b − a) · ωf(d(∆)), f ∈ C([a, b]). (8) Ocena bł ˛edu. Podsumowuj ˛ac, z (6) i (8) otrzymujemy

|I(f) − σ(f, ∆, X∆)| ≤ c · (d(∆))^α, (9) gdzie

- c jest stał ˛a niezale˙zn ˛a od n;

- α = 1, gdy f jest funkcj ˛a monotoniczn ˛a lub f ∈ C¹([a, b]);

- α ∈ (0, 1], gdy f spełnia warunek Höldera z wykładnikiem α.

Mo˙zemy otrzyma´c lepsze oszacowanie, gdy wybierzemy ξi = (xi+ xi−1)/2 (reguła punktu ´srod- kowego).

Przykład. Rozwa˙zamy funkcj ˛e f(x) = x² na przedziale [0, 1]. Oszacuj bł ˛ad (3) w zale˙zno´sci od d(∆), gdy ξ_i= (xi+ xi−1)/2 (i = 1, 2, ..., n).

Wiemy, ˙ze I(f) = 1/3. Zatem I(f ) − σ(f, ∆, X∆) = 1

3−

n i=1

xi+ xi−1

2

(xi− xi−1) =

=

n i=1

1 3

x³_i − x³_i−1

−1

4(xi+ xi−1)²(xi− xi−1)

= 1 12

n i=1

(xi− xi−1)³,

a st ˛ad

|I(f) − σ(f, ∆, X∆)| ≤ 1

12(d(∆))².

(6)

1.2 Przybli˙zone metody - postawienie zadania

W praktyce mamy na ogół mo˙zliwo´s´c obliczania warto´sci funkcji podcałkowej f w pewnych punktach z przedziału [a, b]. Rozwa˙za´c b ˛edziemy metody przybli˙zonego całkowania, w których przybli˙zenia całki I(f) przedstawia si ˛e w postaci kombinacji liniowej takich warto´sci, tj. metody okre´slone wzorami

S(f) =

n i=0

Aif(xi), xi∈ [a, b]. (10)

Wzory te maj ˛a 2n + 3 parametry: liczba składników - n, w ˛ezły - x0, x1, ..., xn i współczynniki - A0, A1, ..., An. W ogólnym przypadku w ˛ezły i współczynniki zale˙z ˛a od n; b ˛edziemy to zaznacza´c pisz ˛ac (o ile to b ˛edzie konieczne) x⁽ⁿ⁾_i i A⁽ⁿ⁾_i (0 ≤ i ≤ n).

Całka I(f ) jest liczb ˛a, a w oparciu o wzory (10) obliczamy przybli˙zenie tej liczby. Zakładamy zbie˙zno´s´c

n i=0

A⁽ⁿ⁾_i f(x⁽ⁿ⁾_i ) →

n→∞I(f).

Oznacza to, ˙ze pomocy tych wzorów mo˙zna obliczy´c I(f) z dowoln ˛a dokładno´sci ˛a, je´sli tylko we´zniemy wystarczaj ˛aco du˙zo w ˛ezłów i współczynników. Z praktycznego punktu widzenia wa˙zna jest szybko´s´c tej zbie˙zno´sci.

W niniejszym opracowaniu, w odniesieniu do wzorów (10), nie b ˛edziemy poruszać problemów o charakterze ogólnym, np. jakie nale˙zy nało˙zyć warunki na w ˛ezły i współczynniki, przy których mamy zapewnion ˛a zbie˙zno´sć. Zwykle zakłada si ˛e, ˙ze f nale˙zy do okre´slonej klasy funkcji, np.

f ∈ C([a, b]), f ∈ C^k([a, b]) (k ≥ 1).

Poni˙zej omówimy dwie proste metody przybli˙zonego całkowania: metod ˛e prostok ˛atów i metod ˛e trapezów.

Zauwa˙zmy, ˙ze suma całkowa σ(f, ∆, X∆) ma posta´c (10): w ˛ezłami s ˛a punkty po´srednie ξ₁, ξ₂, ..., ξn, a współczynnikami odpowiednio ró˙znice x1− x0, x2− x1, ...., xn− xn−1.

1.3 Metoda prostok ˛atów i metoda trapezów Metoda prostok ˛atów. Wzór przybli˙zony

S(f) = b− a n

n i=1

f

a+

i− 1

2

b− a n

nazywamy wzorem prostok ˛atów.

Przedział [a, b] dzielimy punktami xi= a + i(b − a)/n (i = 0, 1, ..., n) na n podprzedziałów o tej samej długo´sci (b−a)/n. Zauwa˙zmy, ˙ze S(f) jest sum ˛a całkow ˛a funkcji f dla ∆ = (x0, x1, ..., xn) i punktów po´srednich ξ₁, ξ₂, ..., ξ_n okre´slonych w nast ˛epuj ˛acy sposób:

ξ_i= xi−1+ xi

2 = a +

i−1

2

b− a n .

Interpretacj ˛e geometryczn ˛a wzoru prostok ˛atów dla n = 3 podano na rysunku.

(7)

Uwaga. Oczywi´scie d(∆) = (b − a)/n. Je˙zeli f ∈ C¹([a, b]), to z (8) wynika nastepuj ˛ace oszacowanie bł ˛edu metody prostok ˛atów

|I(f ) − S(f )| ≤ (b − a)² n M1.

Widzimy, ˙ze je˙zeli f ∈ C¹([a, b]), to powy˙zsze oszacowanie bł ˛edu tej metody maleje ze wzrostem njak 1/n. Poka˙zemy, ˙ze bł ˛ad ten maleje jak 1/n², gdy f ∈ C²([a, b]).

Twierdzenie 1. Je˙zeli f ∈ C²([a, b]), to

I(f) − S(f ) = (b − a)³

24n² f⁽²⁾(η), (11)

gdzie η jest pewnym punktem z przedziału [a, b].

Dowód. Niech x ∈ [xk−1, xk]. Na podstawie wzoru Taylora mamy f(x) = f(ξ_k) + (x − ξ_k)f^′(ξ_k) +1

2(x − ξ_k)²f⁽²⁾(ζ_k),

gdzie ζ_k= ζ_k(x) jest pewn ˛a liczb ˛a le˙z ˛ac ˛a mi ˛edzy x i ξk, na ogół zale˙zn ˛a od x. Całkuj ˛ac powy˙zsz ˛a równo´s´c od x_k−1 do xk otrzymamy

xk

x_k−1

f(x)dx = b− a

n f(ξ_k) + f^′(ξ_k)

xk

x_k−1

(x − ξ_k)dx +1 2

xk

x_k−1

f⁽²⁾(ζ_k(x))(x − ξ_k)²dx. (12)

Drugi składnik po prawej stronie odpada, poniewa˙z

xk

x_k−1

(x − ξ_k)dx = 0. (13)

(8)

Czynnik (x−ξk)²wyra˙zenia podcałkowego trzeciego składnika nie zmienia znaku. Wykorzystuj ˛ac uogólnione twierdzenie o warto´sci ´sredniej dla całek zapiszemy

xk

x_k−1

f⁽²⁾(ζ_k(x))(x − ξ_k)²dx= f⁽²⁾(η_k)

xk

x_k−1

(x − ξ_k)²dx= 1 12

(b − a)³

n³ f⁽²⁾(η_k), (14) gdzie η_k ∈ (xk−1, xk). Uwzgl ˛edniaj ˛ac (12), (13) i (14), dostajemy

I(f) =

n k=1

xk

x_k−1

f(x)dx = b− a n

n k=1

f(ξ_k) + 1 24

(b − a)³ n³

n k=1

f⁽²⁾(η_k) =

= S(f) + 1 24

(b − a)³ n² · 1

n

n k=1

f⁽²⁾(ηk).

Liczba 1 n

n k=1

f⁽²⁾(η_k) jest ´sredni ˛a arytmetyczn ˛a warto´sci drugiej pochodnej w punktach η_k(k = 1, 2, ..., n). Warto´s´c ta znajduje si ˛e mi ˛edzy najmniejsz ˛a i najwi ˛eksz ˛a warto´sci ˛a pochodnej f⁽²⁾(x) na przedziale [a, b], a poniewa˙z pochodna ta jest ci ˛agła na [a, b], wi ˛ec istnieje punkt η ∈ [a, b]

taki, ˙ze

1 n

n k=1

f⁽²⁾(ηk) = f⁽²⁾(η).

Zatem bł ˛ad wzoru prostok ˛atów mo˙zna przedstawi´c w postaci (11). Przykład. Niech f (x) = 1/(1 + x). Przybli˙zona warto´s´c całki I(f) =

1 0

f(x)dx obliczona za pomoc ˛a wzoru prostok ˛atów dla n = 3 wynosi

S(f ) = 1 3

f(1

6) + f (3

6) + f(5 6)

= 478

693 ≈ 0.6897...

Ocenimy teraz dokładno´s´c tego przybli˙zenia. Poniewa˙z f⁽²⁾(x) = 2/(1 + x)³, wi ˛ec istnieje η∈ [0, 1] takie, ˙ze

I(f) −478 693 = 1

108· 1 (1 + η)³. St ˛ad

1

864 ≤ I(f) −478 693 ≤ 1

108. Dokładna warto´s´c całki wynosi ln 2 ≈ 0.69314...

(9)

Metoda trapezów. Wzór przybli˙zony S(f) = b− a

n

f(a) + f (b)

2 +

n−1

i=1

f(a + ib− a n )

nazywamy wzorem trapezów.

Interpretacj ˛e geometryczn ˛a wzoru trapezów dla n = 3 podano na rysunku.

Przedział [a, b] dzielimy punktami xi = a + i(b − a)/n (i = 0, 1, ..., n) na n podprzedziałów o tej samej długo´sci h = (b − a)/n. Pole obszaru domkni ˛etego ograniczonego wykresem funkcji y= f (x) ≥ 0, osi ˛a Ox i prostymi x = a, x = b przybli˙zamy sum ˛a pól n trapezów

S(f ) = h

2(f(x0) + f(x1)) +h

2(f(x1) + f(x2)) + ... +h

2(f (xn−1) + f (xn)) . Twierdzenie 2. Je˙zeli f ∈ C²([a, b]), to

I(f) − S(f) =−(b − a)³

12n² f⁽²⁾(η), (15)

gdzie η jest pewnym punktem z przedziału [a, b].

Dowód. Niech L^k₁ b ˛edzie wielomianem interpolacyjnym funkcji f opartym na punktach xk−1 i xk (patrz rysunek).

(10)

Mamy wówczas

f(x) = L^k₁(x) +f⁽²⁾(ζ_k)

2 (x − xk−1)(x − xk), gdzie ζ_k= ζ_k(x) ∈ (x_k−1, xk). Uwzgl ˛edniaj ˛ac powy˙zsz ˛a równo´s´c, dostajemy

I(f) =

n k=1

xk

x_k−1

f(x)dx =

=

n k=1

xk

x_k−1

L^k₁(x)dx +1 2

n k=1

xk

x_k−1

f⁽²⁾(ζ_k(x))(x − xk−1)(x − xk)dx =

= S(f) + 1 2

n k=1

xk

x_k−1

f⁽²⁾(ζ_k(x))(x − xk−1)(x − xk)dx.

Czynnik (x − x_k−1)(x − xk) w wyra˙zeniach podcałkowych nie zmienia znaku na przedziale [xk−1, xk], mo˙zemy zatem zastosowa´c do powy˙zszych całek uogólnione twierdzenie o warto´sci

´sredniej dla całek. Bior ˛ac pod uwag ˛e, ˙ze

xk

x_k−1

(x − xk−1)(x − xk)dx = −1

6(xk− xk−1)³ mamy

I(f ) = S(f ) +1 2

n k=1

f⁽²⁾(η_k)

xk

x_k−1

(x − xk−1)(x − xk)dx =

= S(f ) −(b − a)³ 12n³

n k=1

f⁽²⁾(η_k) = S(f) − (b − a)³ 12n² · 1

n

n k=1

f⁽²⁾(η_k).

gdzie η_k ∈ (xk−1, xk). Rozumuj ˛ac podobnie jak w metodzie prostok ˛atów stwierdzamy, ˙ze istnieje punkt η ∈ [a, b] taki, ˙ze

1 n

n k=1

f⁽²⁾(η_k) = f⁽²⁾(η).

Zatem bł ˛ad wzoru trapezów mo˙zna przedstawi´c w postaci (15). Podsumowuj ˛ac, metoda prostok ˛atów i metoda trapezów polega na

• równomiernym podziale przedziału całkowania [a, b] na pewn ˛a liczb ˛e podprzedziałów;

• zast ˛apieniu funkcji podcałkowej w ka˙zdym takim podprzedziale wielomianem interpolacyjnym - odpowiedni ˛a stał ˛a lub odpowiednim wielomianem stopnia co najwy˙zej pierwszego.

(11)

Bł ˛ad tych metod maleje ze wzrostem n jak 1/n², gdy f ∈ C²([a, b]). W ka˙zdym podprzedziale mo˙zna stosowa´c interpolacj ˛e wielomianem wy˙zszego stopnia. Na przykład, stosuj ˛ac interpolacj ˛e paraboliczn ˛a mo˙zna otrzyma´c wzory przybli˙zone, dla których bł ˛ad maleje jak 1/n⁴, gdy f ∈ C⁴([a, b]).

1.4 Przypomnienie wybranych faktów

Bł ˛ad interpolacji. Niech f : [a, b] → R i niech Ln b ˛edzie wielomianem interpolacyjnym Lagrange’a funkcji f opartym na w ˛ezłach x0, x1, ..., xn. Je˙zeli f ∈ Cⁿ⁺¹([a, b]) i x ∈ [a, b], to istnieje ξ = ξ_x ∈ (a, b), dla którego

f(x) − Ln(x) = f⁽ⁿ⁺¹⁾(ξ)

(n + 1)! (x − x0)(x − x1)...(x − xn).

Uogólnione twierdzenie o warto´sci ´sredniej dla całek. Je˙zeli na przedziale [a, b] funkcja f jest ci ˛agła, a funkcja g całkowalna i nie zmienia znaku na tym przedziale (tzn. g ≥ 0 lub g ≤ 0), to istnieje punkt ξ ∈ (a, b) taki, ˙ze

b a

f(x)g(x)dx = f (ξ)

b a

g(x)dx.

Moduł ci ˛agło´sci. Niech f : [a, b] → R . Funkcj ˛e ωf : [0, b − a] → R⁺∪ {0}

ωf(δ) = sup

|x−y|≤δ

|f(x) − f(y)|

nazywamy modułem ci ˛agło´sci funkcji f. Własno´sci funkcji ωf :

• lim

δ→0ωf(δ) = 0 ⇔ f jest ci ˛agła na [a, b];

• ωf(δ) ≤ Lδ^α dla f spełniaj ˛acych warunek Höldera ze stał ˛a L i z wykładnikiem α;

• ωf(δ) ≤ M1δ dla f ∈ C¹([a, b]), gdzie M1= max

a≤x≤b|f^′(x)| ;

(12)

2 Zastosowanie wzoru Taylora do oblicze´n przybli˙zonych

Wzór Taylora mo˙zna wykorzysta´c do znajdowania wzorów przybli˙zonych prostej postaci (przybli˙zeniami s ˛a wielomiany), a tak˙ze do przybli˙zonego obliczania warto´sci funkcji. Dla uproszczenia ograniczymy si ˛e do funkcji dwóch zmiennych.

Niech funkcja f (dwóch zmiennych) ma ci ˛agłe pochodne do rz ˛edu (n + 1) wł ˛acznie w pewnym otoczeniu punktu (x0, y0) i niech (x, y) b ˛edzie dowolnym punktem z tego otoczenia. Bierzemy pod uwag ˛e wzór Taylora

f(x, y) = f (x0, y0) + df(x0, y0) +1

2d²f(x0, y0) + ... + 1

n!dⁿf(x0, y0)+

+ 1

(n + 1)!dⁿ⁺¹f(x0+ θ(x − x0), y0+ θ(y − y0)), 0 < θ < 1,

(16)

gdzie d^rf(x^′, y^′) jest ró˙zniczk ˛a r-tego rz ˛edu funkcji f w punkcie (x^′, y^′) liczon ˛a dla przyrostów x− x0 i y − y0:

d^rf(x^′, y^′) =

r k=0

r k

∂^rf

∂x^r−k∂y^k(x^′, y^′) · (x − x0)^r−k(y − y0)^k. W przypadku, gdy x0 = 0 i y0= 0, wzór Taylora nazywamy wzorem Maclaurina.

Wzór Taylora mo˙zna wykorzysta´c do otrzymania wzorów przybli˙zonych dla funkcji f i ocen bł ˛edów, jakie popełniamy stosuj ˛ac te wzory. Po pomini ˛eciu reszty, otrzymujemy wzór przybli˙zony

f(x, y) ≈ f(x0, y0) + df(x0, y0) +1

2d²f(x0, y0) + ... + 1

n!dⁿf(x0, y0). (17) Przybli˙zeniami s ˛a wielomiany zmiennych x i y - punkt (x0, y0) jest ustalony. Dokładniej, d^rf(x0, y0) jest kombinacj ˛a liniow ˛a wielomianów zmiennych x, y postaci (x − x0)^r−k(y − y0)^k (k = 0, 1, ..., r).

Korzystaj ˛ac ze wzoru przybli˙zonego (17) popełniamy bł ˛ad

∆ =

f(x, y) −

f(x0, y0) + df(x0, y0) +1

2d²f(x0, y0) + ... + 1

n!dⁿf(x0, y0)

= 1

(n + 1)!

dⁿ⁺¹f(x0+ θ(x − x0), y0+ θ(y − y0)) ,

który zale˙zy od nieznanego θ ∈ (0, 1). Bł ˛ad ten mo˙zna oszacowa´c, niezale˙znie od θ, nast ˛epuj ˛aco

∆ ≤ Mn+1

(n + 1)!

n+1

k=0

n+ 1 k

|x − x0|^n+1−k|y − y0|^k= Mn+1

(n + 1)!(|x − x0| + |y − y0|)ⁿ⁺¹,

(13)

gdzie Mn+1 spełnia oszacowanie

0≤k≤n+1max sup

θ

∂ⁿ⁺¹f

∂x^n+1−k∂y^k(x0+ θ(x − x0), y0+ θ(y − y0))

≤ Mⁿ⁺¹.

W pobli˙zu punktu (x0, y0) przybli˙zenia te obarczone s ˛a małymi bł ˛edami, które rosn ˛a przy oddalaniu si ˛e od tego punktu.

Wyznaczanie wzorów przybli˙zonych (17) nie zawsze jest łatwe. Je˙zeli funkcja f jest rozwijalna w szereg Taylora w otoczeniu punktu (x0, y0), to mo˙zemy skorzystać z faktu, ˙ze rozwini ˛ecie to jest jednoznaczne. Rozwini ˛ecie to mo˙zemy wyznaczyć jak ˛akolwiek b ˛ad´z metod ˛a, na przykład mo˙zemy wykorzystać tutaj rozwini ˛ecia dla funkcji jednej zmiennej.

Dla funkcji f(x, y) = exp(x + y) mamy rozwini ˛ecie exp(x + y) =

∞ k=0

(x + y)^k

k! dla dowolnych x, y ∈ R. (18)

Mamy tutaj d^kf(0, 0) = (x + y)^k - jest to kombinacja liniowa jednomianów k-tego stopnia zmiennych x i y. Rozwini ˛ecie (18) jest rozwini ˛eciem Maclaurina funkcji f . Bior ˛ac kolejne sumy cz ˛e´sciowe tego rozwini ˛ecia, otrzymujemy wzory przybli˙zone postaci (17)

exp(x + y) ≈ 1, exp(x + y) ≈ 1 + (x + y), exp(x + y) ≈ 1 + (x + y) + 1

2(x + y)², ....

Przykład. Niech s2(x, y) = 1 −1 2

x²+ y²

i s4(x, y) = 1 −1 2

x²+ y² + 1

24(x⁴+ 6x²y²+ y⁴).

Uzasadni´c wzory przybli˙zone

cos(x) cos(y) ≈ s2(x, y), cos(x) cos(y) ≈ s4(x, y), (19) i oceni´c bł ˛ad, jaki popełniamy stosuj ˛ac te wzory.

Przyjmujemy

f(x, y) = cos(x) cos(y).

Skorzystamy z rozwini ˛e´c dla funkcji jednej zmiennej f(x, y) = 1

2(cos(x − y) + cos(x + y)) = 1 2

_∞

k=0

(−1)^k(x − y)^2k (2k)! +

∞ k=0

(−1)^k(x + y)^2k (2k)!

.

Składnikami tych rozwini ˛e´c s ˛a wielomiany zmiennych x i y, i po uporz ˛adkowaniu prawej strony f(x, y) = 1

2

∞ k=0

(−1)^k(x − y)^2k+ (x + y)^2k

(2k)! =

= 1 − 1 2

x²+ y² + 1

24

x⁴+ 6x²y²+ y⁴

− ....

(14)

otrzymamy rozwini ˛ecie Maclaurina funkcji f. Bior ˛ac sumy cz ˛e´sciowe tego rozwini ˛ecia, sum ˛e dwóch pierwszych składników i sum ˛e trzech pierwszych składników, otrzymujemy wzory przybli˙zone (19). Mo˙zemy przyj ˛a´c Mn+1 = 1 dla wszystkich n ≥ 0. Poniewa˙z d³f(0, 0) = 0 i d⁵f(0, 0) = 0, wi ˛ec

|cos(x) cos(y) − si(x, y)| ≤ (|x| + |y|)ⁱ⁺²

(i + 2)! dla i = 2, 4.

W pobli˙zu pocz ˛atku układu współrz ˛ednych przybli˙zenia te obarczone s ˛a małymi bł ˛edami, które rosn ˛a przy oddalaniu si ˛e od tego punktu (patrz rysunek).

Do obliczania przybli˙zonych warto´sci skomplikowanych wyra˙ze´n cz ˛esto stosuje si ˛e wzory przybli˙zone oparte na pierwszej ró˙zniczce

f(x, y) ≈ f(x0, y0) + df (x0, y0) = f (x0, y0) +∂f

∂x(x0, y0) · (x − x0) +∂f

∂y(x0, y0) · (y − y0). (20) Przykład. Korzystaj ˛ac ze wzoru przybli˙zonego (20), obliczy´c przybli˙zon ˛a warto´s´c wyra˙zenia

(2.02)³+ (0.97)⁴.

Przyjmujemy tutaj f(x, y) =

x³+ y⁴ oraz x0= 2 i y0 = 1. Mamy nast ˛epnie

∂f

∂x = 3x² 2

x³+ y⁴, ∂f

∂y = 2y³ x³+ y⁴.

Korzystaj ˛ac z przybli˙zenia (20), otrzymamy f(2.02, 0.97) =

(2.02)³+ (0.97)⁴ ≈ 3 + 2 · 0.02 −2

3 · 0.03 = 3.02 . Dokładna warto´s´c wyra˙zenia wynosi 3.0212...

(15)

3 Poszukiwanie ekstremów funkcji

Rozwa˙zamy problem wyznaczania ekstremów (lokalnych) funkcji wielu zmiennych. Analityczne metody, omawiane w kursie analizy matematycznej, pozwalaj ˛a na efektywne wyznaczenie ek- stremów jedynie prostych funkcji. Trudno´sci, na które natrafiamy przy wyznaczaniu ekstremów doprowadziły do powstania wielu metod przybli˙zonego wyznaczania ich. Jedn ˛a z najprostszych metod jest metoda najszybszego spadku, któr ˛a obecnie przedstawimy. Zało˙zymy, ˙ze poszukujemy minimum (podobnie poszukuje si ˛e maksimum).

Poprzedzaj ˛ac opis metody najszybszego spadku, przedstawimy istot ˛e jej działania na przykładzie funkcji dwóch zmiennych. Niech funkcja f ma minimum w punkcie M. Rozwa˙zmy warstwic ˛e funkcji f znajduj ˛ac ˛a si ˛e w pewnym otoczeniu punktu M i niech M0b ˛edzie punktem tej warstwicy.

Wektor −→n = − grad f(M0) wskazuje kierunek najszybszego spadku funkcji f w punkcie M0. Zatem −→n mo˙ze by´c uznany jako odpowiedni kierunek przy poszukiwania punktu M, gdy znaj- dujemy si ˛e w punkcie M0.

Metoda najszybszego spadku. Niech f : D → R, gdzie D ⊂ Rⁿ.Zakładamy, ˙ze f ∈ C¹(D).

Opis metody: kolejne przybli˙zenia M0, M1, M2, ....wyznaczamy według schematu (0) Ustal przybli˙zenie pocz ˛atkowe M0.Przyjmij i = 0.

(1) Oblicz warto´s´c funkcji f (Mi) i wyznacz kierunek poszukiwa´n −→ni = −grad f (Mi) w punkcie Mi.

(2) Wzdłu˙z kierunku −→niokre´sl λi>0 minimalizuj ˛ace warto´s´c funkcji f(Mi+λi−→ni) i współrz ˛edne nast ˛epnego punktu

Mi+1 = Mi+ λi−→ni. (3) Podstaw i := i + 1 i powtórz cykl od punktu (1).

Ci ˛ag punktów (Mi) (o ile proces obliczeniowy jest niesko´nczony) mo˙ze nie by´c zbie˙zny do punktu, w którym f ma minimum (przyjmuj ˛ac, ˙ze taki punkt istnieje). Je˙zeli f jest ´sci´sle wypukła w D

(16)

i D jest zbiorem wypukłym, to globalne minimum f wyst ˛epuje tylko w jednym punkcie i ci ˛ag (Mi) jest zbie˙zny do tego punktu.

Stosuje si ˛e ró˙zne warunki zako´nczenia oblicze´n, np. wyznaczony punkt jest punktem stacjonarnym (punktem bliskim takiemu punktowi) funkcji f. Przebieg powy˙zszego algorytmu przedstawiono na poni˙zszym rysunku (na przykładzie funkcji dwóch zmiennych).

Przykład. Rozwa˙zamy funkcj ˛e f(x, y) = x³+y³−6xy +10. Funkcja ta ma minimum w punkcie M(2, 2), przy czym f(M) = 2. Wyznaczymy teraz 1-sze przybli˙zenie M1 punktu M metod ˛a najszybszego spadku, przyjmuj ˛ac przybli˙zenie pocz ˛atkowe M0(1, 2).

Obliczamy

∂f

∂x = 3x²− 6y, ∂f

∂y = 3y²− 6x i grad f(M0) = [−9, 6].

Wektor −→n0 = [9, −6] wskazuje kierunek najszybszego spadku funkcji f w punkcie M0. Wzdłu˙z tego kierunku poszukujemy punktu M1 = M0+ λ0· −→n0 = (1 + 9λ0,2 − 6λ0), przy czym λ0 >0 ma minimalizowa´c warto´s´c funkcji f(1 + 9λ0,2 − 6λ0)

f(1 + 9λ0,2 − 6λ0) = 513λ³₀+ 783λ²₀− 117λ0+ 7.

Powy˙zszy wielomian ma minimum w punkcie λ0 ≈ 0.0699.., i st ˛ad ostatecznie M1(1.629.., 1.5805..).

(17)

Otrzymali´smy f(M1) ≈ 2.822..., podczas gdy f(M0) = 7. Powtarzaj ˛ac obliczenia b ˛edziemy otrzymywa´c coraz lepsze przybli˙zenia punktu M.

Kierunek −→n0 oraz punkty M0 i M1 zostały zaznaczone na zał ˛aczonym planie warstwicowym funkcji f .

Metoda najszybszego spadku jest jedn ˛a z najprostszych metod przybli˙zonego wyznaczania ek- stremów funkcji. Nie nale˙zy ona do metod ”najefektywniejszych” - wad ˛a tej metody jest wolna zbie˙zno´s´c kolejnych przybli˙ze´n M0, M1, M2, ....

Metoda najszybszego spadku jest metod ˛a gradientow ˛a, metod ˛a, w której korzystamy zarówno z informacji o warto´sci funkcji jak i warto´sci jej gradientu. Opracowano ”efektywne” metody gradientowe, du˙zo bardziej ”zło˙zone” ni˙z metoda najszybszego spadku.

(18)

4 Metoda najmniejszych kwadratów

Załó˙zmy, ˙ze aby wyznaczy´c pewn ˛a zale˙zno´s´c fizyczn ˛a y= f(x)

dokonali´smy pomiarów wielko´sci x i y. Na podstawie tych pomiarów, cz ˛esto obarczonych bł ˛e- dami, staramy si ˛e odtworzy´c f - oczywi´scie w przybli˙zonej postaci. B ˛edziemy poszukiwa´c funkcji F, funkcji aproksymuj ˛acej (przybli˙zaj ˛acej) f, która ”wygładzi” bł ˛edy pomiarowe.

Danych jest (n+1) ró˙znych punktów x0, x1, ..., xnz przedziału [a, b] oraz warto´sci pewnej funkcji y= f (x) w tych punktach

y0= f(x0), y1 = f(x1), ..., yn= f(xn).

Przyjmujemy, ˙ze funkcja przybli˙zaj ˛aca F jest kombinacj ˛a liniow ˛a wybranych funkcji z niez- nanymi współczynnikami:

F(x) =

m j=0

ajϕ_j(x), x ∈ [a, b], (21)

gdzie ϕ₀, ϕ₁, ..., ϕm s ˛a funkcjami okre´slonymi na przedziale [a, b]. Funkcje te s ˛a dane, zadane na podstawie wiedzy o badanym zjawisku.

W metodzie najmniejszych kwadratów, współczynniki a0, a1, ..., am dobieramy tak, ˙zeby funkcja H: R^m+1 → R⁺∪ {0},

H(a0, a1, ..., am) =

n i=0

(yi− F (xi))² =

n i=0



yi−

m j=0

ajϕj(xi)





2

, (22)

osi ˛agała warto´s´c minimaln ˛a.

Funkcj ˛e F , dla której H osi ˛aga warto´s´c minimaln ˛a nazywamy funkcj ˛a optymaln ˛a (wzgl ˛edem układu (ϕ₀, ϕ₁, ..., ϕ_m) na zbiorze {x0, x1, ..., xn}).

Istota metody najmniejszych kwadratów. Zakładamy, ˙ze funkcja przybli˙zaj ˛aca F ma posta´c (21). Licz ˛ac si ˛e z tym, ˙ze bł ˛edy

δi = |yi− F (xi)| (i = 0, 1, ..., n)

s ˛a nieuniknione przy ka˙zdym wyborze współczynników a0, a1, ..., am, współczynniki te dobieramy tak, aby suma kwadratów bł ˛edów

δ²₀+ δ²₁+ ... + δ²_n

była mo˙zliwie mała (st ˛ad nazwa metody). Innymi słowy, za najbardziej zgodne z wynikami pomiaru uwa˙zamy te warto´sci a0, a1, ..., am, dla których H = H(a0, a1, ..., am) ma najmniejsz ˛a warto´s´c.

(19)

Aby znale´z´c współczynniki a0, a1, ..., amobliczamy pochodne cz ˛astkowe ∂H/∂ak(k = 0, 1, ..., m).

Z warunku zerowania si ˛e tych pochodnych (patrz metody poszukiwania ekstremów funkcji wielu zmiennych):

∂H

∂ak = −2

n i=0



yi−

m j=0

ajϕ_j(xi)



 ϕ_k(xi) = 0, k= 0, 1, ..., m,

otrzymamy układ (m + 1) równa´n liniowych z (m + 1) niewiadomymi a0, a1, ..., am

m j=0

_n

i=0

ϕj(xi)ϕk(xi)

· aj =

n i=0

ϕk(xi)yi, k= 0, 1, ..., m, (23) zwany układem równa´n normalnych.

Twierdzenie. Je˙zeli macierz układu (23) jest nieosobliwa, to istnieje dokładnie jedna funkcja optymalna F postaci (21). Rozwi ˛azanie a0, a1, ..., am układu (23) jest ci ˛agiem współczynników tej funkcji.

Dowód tego twierdzenia pominiemy.

Trzy szczególne przypadki.

• Funkcja przybli˙zaj ˛aca F jest stał ˛a

F(x) = a0. St ˛ad

H(a0) =

n i=0

(yi− a0)² i H^′(a0) = −2

n i=0

(yi− a0) = 0.

Optymalna stała

a0 = 1 n+ 1

n i=0

yi (24)

jest ´sredni ˛a arytmetyczn ˛a (n + 1) warto´sci funkcji f.

(20)

• Funkcja przybli˙zaj ˛aca F jest postaci

F(x) = a0x.

St ˛ad

H(a0) =

n i=0

(yi− a0xi)² i H^′(a0) = −2

n i=0

(yi− a0xi) · xi= 0.

F jest optymaln ˛a funkcj ˛a, gdy

a0 =

n i=0

xiyi

n i=0

x²_i

. (25)

• Funkcja przybli˙zaj ˛aca F jest postaci

F(x) = a0+ a1x.

St ˛ad

H(a0, a1) =

n i=0

(yi− (a0+ a1xi))²

oraz

∂H

∂a0 = −2

n i=0

(yi− (a0+ a1xi)) = 0

∂H

∂a1 = −2

n i=0

(yi− (a0+ a1xi)) · xi = 0

Układ równa´n normalnych ma posta´c











(n + 1) · a0+ _n

i=0

xi

· a1 =

n i=0

yi

_n

i=0

xi

· a0+ _n

i=0

x²_i

· a1 =

n i=0

yixi

.

Macierz tego układu jest nieosobliwa, poniewa˙z

(n + 1)

n i=0

x²_i − _n

i=0

xi

2

=

n i,j=0(i<j)

(xi− xj)² = 0.

(21)

Przykład. W wyniku osiowego rozci ˛agania próbki stalowej w maszynie wytrzymało´sciowej otrzymano dane eksperymentalne (ε, σ) (Glinicka A., IL PW) :

(0, 0), (0.01, 19.6), (0.024, 39.2), (0.033, 58.8), (0.04, 78.5), (0.053, 98), (0.066, 117.7), (0.071, 127.5), (0.079, 137.3),

gdzie σ - napr ˛e˙zenie [MPa] oraz ε - wydłu˙zenie wzgl ˛edne [%]. Na podstawie tych 9 pomiarów wyznaczy´c zale˙zno´s´c mi ˛edzy σ i ε.

Wiemy, ˙ze dla małych odkształceń i napr ˛e˙zeń, na mocy prawa Hooke’a, zwi ˛azek mi ˛edzy nimi jest liniowy. Taki charakter maj ˛a powy˙zsze dane. Z tego powodu do aproksymacji nale˙zy wybrać funkcj ˛e liniow ˛a postaci F (ε) = a · ε, gdzie a jest nieznanym parametrem. Współczynnik proporcjonalno´sci a zwany jest modułem Younga.

F jest optymaln ˛a funkcj ˛a, gdy współczynnik a ma posta´c (25). Dla powy˙zszych danych ekspery- mentalnych otrzymujemy a = 1.792 · 10⁵ MPa. Dane pomiarowe (ε, σ) oraz wykres optymalnej funkcji przedstawiono na powy˙zszym rysunku.

Metoda najmniejszych kwadratów jest najpopularniejsz ˛a metod ˛a opracowania danych pomi- arowych. Metoda ta została wprowadzona przez A.M.Lagrange’a w 1805. W latach 1809-1810 K.F.Gauss, który twierdził, ˙ze u˙zywał jej od 1794 r., przedstawił podstawy rachunkowe tej metody i podał dla niej uzasadnienie probabilistyczne.