Metody numeryczne- wybrane zagadnienia

dr Bogdan Grabiec

Zielona Góra 2010

Tytuł projektu:

Unowocześnienie programów, poprawa jakości kształcenia oraz otwarcie nowej
specjalności ekofizyka na kierunku fizyka w Uniwersytecie Zielonogórskim,

Wniosek:

POKL 04.01.01-00-041/09-00
Program Operacyjny Kapitał Ludzki

Priorytet:

IV. Szkolnictwo wyższe i nauka

Działanie:

4.1 Wzmocnienie i rozwój potencjału dydaktycznego uczelni oraz zwiększenie liczby
absolwentów kierunków o kluczowym znaczeniu dla gospodarki opartej na wiedzy,

Poddziałanie:

4.1.1 Wzmocnienie potencjału dydaktycznego uczelni

Beneficjent:

Uniwersytet Zielonogórski, ul. Licealna 9, 65-417 Zielona Góra

Spis treści

1 Interpolacja. Interpolacja Lagrange’a

3

2 Aproksymacja

9

2.0.1

Aproksymacja średniokwadratowa . . . . . . . . . . . .

9

2.0.2

Aproksymacja jednostajna . . . . . . . . . . . . . . . . 11

3 Numeryczne różniczkowanie

15

4 Całkowanie numeryczne

17

4.0.3

Metoda trapezów . . . . . . . . . . . . . . . . . . . . . 18

4.0.4

Metoda Simpsona . . . . . . . . . . . . . . . . . . . . . 20

5 Numeryczne rozwiązywanie równań

23

6 Równania rożniczkowe

27

1

2

SPIS TREŚCI

Rozdział 1

Interpolacja. Interpolacja
Lagrange’a

Wyobraźmy sobie że mamy zadaną funkcje y=f(x) dla x

0

∈ [x

1

, x

2

]. Określe-

nie wartości funkcji dla x z tego przedziału, nie stanowi problemu, wystarczy
wyznaczyć y

0

= f (x

0

). Sprawa się komplikuje gdy nie dysponujemy już ścisłą

postacią funkcji a jedynie zbiorem punktów (x

1

, y

1

), ..., (x

N

, y

n

). Sposób wy-

znaczenia przybliżonej wartości funkcji dla argumentu x

0

znajdującego sie

między zadanymi punktami już nie jest takie proste. Jednym z sposobów
realizujących powyższe zagdanienie nosi nazwę interpolacji.

Przykład 1

3

4

ROZDZIAŁ 1. INTERPOLACJA. INTERPOLACJA LAGRANGE’A

Dysponując dwoma punktami o współrzędnych : (x

1

, y

1

= f (x

1

)) oraz

(x

2

, y

2

= f (x

2

)).Określimy przybliżoną watość funkcji dla argumentu x

0

z

przedziału należącego do [x

1

, x

2

]. Problem ten możemy rozwiązać następują-

co: szukamy równanie prostej przechodzącej przez dane punkty

y

1

− y

x

1

− x

=

y

2

− y

1

x

2

− x

1

.

Na podstawie powyższego równania, możemy wyznaczyć wartość y

y =

y

2

− y

1

x

2

− x

1

(x − x

1

) + y

1

a tym samym, przybliżona wartość dla x=x0

y

0

=

y

2

− y

1

x

2

− x

1

(x

0

− x

1

) + y

1

Przykład 2

Zagadnienie to możemy rozszerzyć o znajomość wspólrzędnych następnego
punktu. Tym razem zadamy sobie pytanie jaki wielomian, tzn. o jakim stop-
niu możemy przeprowadzi˙c przez podane trzy punkty. Odpowiedzią na zada-
ny problem jest funkacja kwadratowa postaci f (x) = a

2

x

2

+a

1

x+a

0

. Musimy

znaleźć wielkości a,b i c wyrażone poprzez znajomość x

1

, x

2

, x

3

, y

1

, y

2

, y

3

Uzy-

skujemy układ równań postaci:











a

2

x

2

1

+ a

1

x

1

+ a

0

= y

1

a

2

x

2

2

+ a

1

x

2

+ a

0

= y

2

a

2

x

2

3

+ a

1

x

3

+ a

0

= y

3

Rozwiązanie ma postać:

a2 =

x

2

1

x

1

y

1

x

2

2

x

2

y

2

x

2

3

x

3

y

3

W

a1 =

x

2

1

y

1

1

x

2

2

y

2

1

x

2

3

y

3

1

W

a2 =

y

1

x

1

1

y

2

x

2

1

y

3

x

3

1

W

5

gdzie

W =

x

2

1

x

1

1

x

2

2

x

2

1

x

2

3

x

3

1

= (x

3

− x

2

)(x

3

− x

1

)(x

2

− x

1

)

Co pozwala zapisać ostatecznie:

W

3

(x) = y

1

(x − x

2

)(x − x

3

)

(x2 − x1)(x3 − x1)

−y

2

(x − x

1

)(x − x

2

)

(x

2

− x

1

)(x

3

− x

2

)

+y

3

(x − x

1

)(x − x

2

)

(x

3

− x

1

)(x

3

− x

1

)

Możemy sformułować następującą własność w postaci twierdzenia.

Twierdzenie 1.0.1 Dla zadanych N punktów :(x

1

, y

1

), ...., (x

N

, y

N

), które

będziemy nazywać węzłami interpolacyjnymi, istnieje dokładnie jeden wie-
lomian stopnia N, postaci

W

N

(x) =

N

X

i

=0

a

i

x

i

.

dowód
Zadane punkty muszą spełniać powyższe r ˙ownania











a

0

+ a

1

x

1

+ ... + a

N

x

N

1

= y

1

...........................
a

0

+ a

1

x

N

+ ... + a

N

x

N

N

= y

N

6

ROZDZIAŁ 1. INTERPOLACJA. INTERPOLACJA LAGRANGE’A

Ponownie, rozwiązania dla powyższego układu stosując twierdzenie Cra-

mera, możemy zapisać w postaci:

a

i

=

1

D

n

X

j

=0

y

j

D

ij

gdzie D

ij

oznacza dopełnienia algebraiczne elementów i-tej kolumny macierzy

A postaci

A =








1

x

1

x

2

1

... x

N

1

1

x

2

x

2

2

... x

N

2

... ...

...

... ...

1

x

N

x

2

N

... x

N

N








co kończy dowód. Warto wspomnieć że wyznacznik macierzy A nie jest trudny
do obliczenia.

|A| =

Y

0¬j<i¬n

(x

i

− x

j

) 6= 0

Jest to tzw. wyznacznik Vandermonde’a.

Twierdzenie 1.0.2 Dla danych węzłów interpolacyjnych x

0

, x

1

, ..., x

N

∈ [a, b]

i danej funkcji

f ∈ C

[n+1]

[a, b] możemy dla dowolnego punktu x ∈ [a, b] zna-

leźć takie

ζ = ζ(x) ∈ [a, b] że

f (x) = P (x) +

f

(n+1)

(ζ(x))

(n + 1)!

(x − x

0

)(x − x

1

)...(x − x

N

)

gdzie

P (x) =

N

X

k

=0

f (x

k

)L

n,k

(x),

L

n,k

(x) =

N

Y

i

=0,i6=k

(x − x

i

)

(x

k

− x

i

)

dla k = 0, 1, ..., n.
dowód:
Jeśli x = x

k

dla k = 0, 1, ..., n, f (x

k

) = P (x

k

).

Natomiast gdy x 6= x

k

wprowadzam funkcję:

g(t) = f (t) − P (t) − [f (x) − P (x)]

(t − x

0

)(t − x

1

)...(t − x

N

)

(x − x

0

)(x − x

1

)...(x − x

N

)

= f (t) − P (t) − [f (x) − P (x)]

n

Y

i

=0

(t − x

i

)

(x − x

i

)

7

Patrząc na funkcję g, widać że musi być klasy c

[n+1]

[a, b]. Dla t = x

k

oraz

t = x g = 0. W związku z tym mogę wykorzystać n+1 razy twierdzenie
Rolle’a. Wobec tego mogę zapisać:

0 = g

(n+1)

(ζ) = f

(n+1)

(ζ) − P

(n+1)

(ζ) − [f (x) − P (x)]

d

(n+1)

dt

(n+1)

n

Y

i

=0

(t − x

i

)

(x − x

i

)

|

t

=ζ

.

Wielkość P (x) jest wielomianem stopnia n. Pochodna (n+1) rzędu z wielkości
daje nam w wyniku zero. Natomiast

d

(n+1)

dt

(n+1)

n

Y

i

=0

(t − x

i

)

(x − x

i

)

|

t

=ζ

=

(n + 1)!

Q

n
i

=0

(x − x

i

)

.

ostatecznie uzyskujemy

0 = f

(n+1)

(ζ) − 0 − [f (x) − P (x)]

(n + 1)!

Q

n
i

=0

(x − x

i

)

.

f (x) = P (x) +

f

(n+1)

(ζ)

(n + 1)!

n

Y

i

=0

(x − x

i

).

Przykład 3

Dla danych punktów:(-2,3),(1,1),(2,-3),(4,8) znaleźć wielomian interpolacyj-
ny
rozwiązanie

W

4

(x) = 3

(x − 1)(x − 2)(x − 4)

(−2 − 1)(−2 − 2)(−2 − 4)

+ 1

(x + 2)(x − 2)(x − 4)

(1 − 3)(1 − 2)(1 − 4)

− 3

(x + 2)(x − 1)(x − 4)

(2 + 2)(2 − 1)(2 − 4)

+ 8

(x + 2)(x − 1)(x − 1)

(4 + 2)(4 − 1)(4 − 2)

.

Oszacowanie błędu wzoru interpolacyjnego. Problem tego zagadnienia po-

lega na wyznaczeniu wielkości

ε(x) = f (x) − W

n

(x),

która będzie nas informowała o tym jak poważne są odchylenia wartości funk-
cji od wartości opartej na wielomianie interpolacyjnym. W tym celu wpro-
wadzamy funkcje typu:

ϕ(u) = f (u) − W

n

(u) + K(U − x

0

)(u − x

1

)...(u − x

n

),

wykorzystujac twierdzenie Rolle’a (n+1) razy, uzyskamy

K =

f

(n+1)

(ζ)

(n + 1)!

a tym samym

|f (x) − W

n

(x)| ¬

M

n

+1

(n + 1)!

|W

n

(x)|.

8

ROZDZIAŁ 1. INTERPOLACJA. INTERPOLACJA LAGRANGE’A

Rozdział 2

Aproksymacja

Wprowadzmy następujące oznaczenia:
f(x) - funkcja która chcemy przbliżyć,
F(x) - funkcja która przybliżamy funkcję f(x) i nazywamy ją funkcją aprok-
symującą.

Warto wspomieć że jeśli aproksymacje rozważamy na przedziale to bę-

dziemy ją nazywać aproksymacją intergralna, natomiast jeśli na zbiorze dys-
kretnym to aproksymacją punktową.

Poprzez aproksymacje funkcji f(x) rozumiemy znalezienie takiej funkcji

F(x)

F (x) = a

0

φ

0

(x) + ... + a

m

φ

m

(x)

gdzie φ

0

, φ

1

, ..., φ

m

oznacza bazę w przestrzeni m+1 wymiarowej, tak aby np.

minimalizowała normę ||f (x) − F (x)||. Całe zagadnienie sprowadza się do
wyznaczenia współczynników a

0

, a

1

, ..., a

m

.

2.0.1

Aproksymacja średniokwadratowa

Dla tego problemu rozważymy aproksymacje punktową, której podstawą jest
zminmalizowanie normy postaci

||F (x) − f (x)|| =

n

X

i

=0

(F (x

i

) − f (x

i

))

2

gdzie i = 1, 2, 3, ...n.

Jeśli rozważymy wykorzystanie bazy przestrzeni w postaci:

φ

0

(x) = 1, φ

1

(x) = x, φ

2

(x) = x

2

, ...., φ

m

(x) = x

m

.

Funkcje aproksymujacą możemy zapisać w postaci

F (x) =

m

X

i

=0

φ

i

(x)x

i

.

9

10

ROZDZIAŁ 2. APROKSYMACJA

Definiujemy funkcje H postaci:

H(a

0

, a

1

, a

2

, ..., a

m

) =

n

X

j

=1

f

j

−

m

X

i

=0

a

i

φ

i

(x

j

)

!

2

Warunek konieczny pozwalający znaleźć minimum:

∂H
∂a

k

= −2

n

X

j

=0

f

j

−

m

X

i

=0

a

i

φ

i

(x

j

)

!

x

k
j

= 0

Co pozwala zapisać powyższa zależność w postaci

m

X

i

=0

a

i

g

ik

= ρ

k

, k = 0, 1, 2, ..., m

gdzie

g

ik

=

n

X

j

=0

x

i

+k

j

ρ

k

=

n

X

j

=0

f (x

j

)x

k
j

Przykład: regresja liniowa

Dysponujemy znajmościa współrzędnych punktów:(x

i

, y

i

), i = 1, 2, ..., m

Aproksymacji dokonujemy za pomocą funkcji:

F (x) = a

0

+ a

1

x.

W tym przypadku musimy zminimalizowac wyrażenie:

H(a

0

, a

1

) =

m

X

j

=1

(f (x

j

) − (a

0

+ a

1

x

j

))

2

.

Warunek konieczny na istnienie minimum:























∂H

(a

0

,a

1

)

∂a

0

= −2

m

X

j

=1

(f (x

j

) − (a

0

+ a

1

x

j

)) = 0

∂H

(a

0

,a

1

)

∂a

1

= −2

m

X

j

=1

x

j

(f (x

j

) − (a

0

+ a

1

x)) = 0

Rozwiązując powyższe równania uzyskujemy

a

0

=

hyihx

2

i − hxihxyi

hx

2

i − (hxi)

2

a

1

=

hxyi − hxihyi

hx

2

i − (hxi)

2

11

gdzie

hxi =

m

X

i

=0

x

i

m

hyi =

m

X

i

=0

y

i

m

hxyi =

m

X

i

=0

x

i

y

i

m

hx

2

i =

m

X

i

=0

(x

i

)

2

m

2.0.2

Aproksymacja jednostajna

W tym przypadku rozważamy funkcję f (x), x ∈< a, b >. Szukamy najmniej-
szej maksymalnej różnicy

||F (x) − f (x)|| = sup

x∈<a,b>

|F (x) − f (x)|

Metoda szeregów potęgowych

Metoda ta jest jedną z najbardziej znanych sposobów aproksymacji jedno-
stajnej. Polega ona na przybliżaniu funkcji szregiem potęgowym, oczywiście
skończonym.

F (x) = f (x

0

) +

f

(1)

(x

0

)

1!

(x − x

0

) + ... +

f

(m)

(x

0

)

m!

(x − x

0

)

m

W tym przypadku błąd aproksymacji wyznaczymy za pomocą zależności

|f (x) − F (x)| ¬

M

m

+1

(m + 1)!

δ

m

+1

Przykład

Funkcje f (x) = sin(x), dla x ∈< −0.1; 0.1 > możemy przybliżyć funkcją

sin(x) = x −

x

3

6

+

x

5

120

12

ROZDZIAŁ 2. APROKSYMACJA

Przybliżenie Pade’go

W tym przypadku funkcja aproksymujacą przyjmuje postać

F (x) = R

n,k

(x) =

L

n

(x)

M

k

(x)

=

a

0

+ a

a

x + ... + a

n

x

n

1 + b

1

x + ... + b

k

x

k

gdzie N = n + m. Na podaną funkcję nakładamy następującę warunki

1. f (0) = R

n,k

(0) =

L

n

(0)

M

k

(0)

2. pochodne f

(m)

(0) − R

(m)
n,k

(0) = 0 dla m=1,2,..., N

Rozwinięcie funkcji f (x) w szereg Taylora w x = 0 można zapisać:

f (x) =

∞

X

i

=0

c

i

x

i

.

Błąd aproksymacji będzie miał następującą postać:

ε(x) = f (x) −

L

n

(x)

M

k

(x)

.

Uwzględniając warunki (1) i (2) uzyskujemy równanie:

∞

X

i

=0

c

i

x

i

k

X

i

=0

b

i

x

i

−

n

X

i

=0

a

i

x

i

=

∞

X

i

=n+k+1

d

i

x

i

które prowadzi do następujących zależności pomiędzy wspołczynnikami a, b
i c















































c

j

= = 0

j < 0

b

j

= = 0

j > k

k

X

j

=0

c

n

+k−s−j

b

j

= 0

s = 0, 1, ..., k − 1

a

r

=

k

X

j

=0

c

r−j

b

j

r = 0, 1, ..., n.

Przykład

Rozważmy funkcje f (x) = e

x

na przedziale x ∈< −0.5; 0.5 >. Znajdziemy

przybliżenie Pade’go R

2,1

(x), n = 2, k = 1, N = 3 dla tej funkcji.

Rozwinięcie Taylora ma postać

f (x) = 1 + x +

x

2

2

+

x

3

6

= R

3,0

(x)

13

wobec tego c

0

= 1, c

1

= 1, c

2

= 0.5, c

3

= 1/6.

Z równań na wspołczynniki b i c otrzymujemy:

c

3−0−0

b

0

+ c

3−0−1

b

1

= c

3

b

0

+ c

2

b

1

=

1
6

+

1
2

b

1

= 0 ⇒ b

1

= −

1
3

a

0

= c

0

b

0

= 1, a

1

= c

1

b

0

+ c

0

b

1

= 1 −

1
3

=

2
3

, a

2

= c

2

b

0

+ c

1

b

1

=

1
2

−

1
3

=

1
6

.

Ostatecznie możemy zapisać

R

2,1

=

1 +

2
3

x +

1
6

x

2

1 −

1
3

x

.

14

ROZDZIAŁ 2. APROKSYMACJA

Rozdział 3

Numeryczne różniczkowanie

Bogata w informację jest znajomość pochodnej pewnej wielkości. W fizyce
najczęściej podawane przykłady to definicje prędkości i przyśpieszenia.

v(t) =

dr

(t)

dt

a(t) =

dv

(t)

dt

=

d

2

r

(t)

dt

2

Znając równanie toru ruchu ciała, czyli funkcje r = r(t), jesteśmy w sta-
nie powiedzieć z jaką szybkością porusza się ciało oraz jaka działa na niego
siła. Pierwszą uzyskamy z wyznaczenia pochodnej polożenia po czasie a dru-
gą z pomnożenia masy przez drugą pochodną polożenia po czasie. Na tym
przykładzie wida˙c, że, umiejętność policzenia pochodnej z danej funkcji jest
niezbędna dla każdego inżyniera, fizyka czy matematyka.

1. Pochodna numeryczna dwu-punktowa

Rozważmy funkcję f ∈ C

2

[a, b], oraz dwa punkty x

0

, x

1

= x

0

+ h ∈ [a, b].

Korzystając z twierdzenia z rodziału 1, możemy zapisać

f (x) = f

0

x − x

1

x

0

− x

1

+ f

1

x − x

0

x

1

− x

0

+

(x − x

0

)(x − x

1

)

2!

f

”

(ζ)

Po podstawieniu i obliczeniu pochodnej po x

f

′

(x) =

f (x

0

+ h) − f (x

0

)

h

+((x−x

0

)−h)f

”

(ζ(x))+

(x − x

0

)(x − x

0

− h)

2!

d

dx

f

”

(ζ(x)).

Gdy podstawimy za x = x

0

, uzyskamy

f

′

(x

0

) =

f (x

0

+ h) − f (x

0

)

h

−

h

2

f

”

(ζ).

czyli

f

′

(x

0

) =

f (x

0

+ h) − f (x

0

)

h

+ O(h).

15

16

ROZDZIAŁ 3. NUMERYCZNE RÓŻNICZKOWANIE

W podobny sposób możemy wyprowadzić wzory:

2. Pochodna numeryczna trzy-punktowa

f

′

(x

0

) =

f (x

0

+ h) − f (x

0

− h)

2h

+ O(h

2

)).

3. Pochodna numeryczna piecio-punktowa

f

′

(x

0

) =

f (x

0

− 2h) − 8f (x

0

− h) + 8f (x

0

+ h) + f (x

0

+ 2h)

12h

+ O(h

4

).

Natomiast wyższe pochodne wyznaczymy za pomocą wzorów:
5. Druga pochodna numeryczna trzy-punktowa

f ”(x

0

) =

f (x

0

+ h) − 2f (x

0

) + f (x

0

− h)

h

2

+ O(h

2

).

oraz
6. Druga pochodna numeryczna piecio-punktowa

f ”(x

0

) =

−f (x

0

− 2h) + 16f (x

0

− h) − 30f (x

0

) + 16f (x

0

+ h) − f (x

0

+ 2h)

12h

2

+O(h

2

).

Powyższe równania można wyprowadzić korzystając z rozwinięcia w sze-

reg Taylora funkcji f (x) w otoczeniu punktu x

0

f (x) = (f (x

0

) +

f

′

(x

0

)

1

(x − x

0

) +

f “(x

0

)

2!

(x − x

0

)

2

+ ... +

f

(n)

(x

0

)

n!

(x − x

0

)

n

+ ..

Dla x = x

k

+1

oraz wiedząc, ze x

k

+1

− x

0

= x

k

+1

− x

k

= h

f (x

k

+1

) = f

k

+1

= f

k

+ f

′

k

· h +

1
2

f

′′

k

· h

2

+

1
6

f

(3)

k

· h

3

+ ...

Przykład
Wyprowadzimy powyższa metoą wzór na pierwszą pochodną numeryczna
3-punktowa. Wymagana jest zanjomość:

f (x

k

+1

) = f

k

+1

= f

k

+ f

′

k

· h +

1
2

f

′′

k

· h

2

+ O(h

3

)

f (x

k−

1

) = f

k−

1

= f

k

− f

′

k

· h +

1
2

f

′′

k

· h

2

+ O(h

3

)

Tworzymy następującą kombinacje:

f

k

+1

− f

k−

1

= f

′

k

· h + O(h

3

)

Po przeksztłceniu uzyskujemy:

f

′

k =

f

k

+1

− f

k−

1

h

+ O(h

2

).

Rozdział 4

Całkowanie numeryczne

Sytuacja jest bardzo podobna jak do analizowanych zagadnień w rodziałe
poprzednim. Tym razem szukamy sposobu na obliczenie

r(t) =

Z

t

2

t

1

v(t)dt

lub

v(t) =

Z

t

2

t

1

a(t)dt.

Zagadnienie możemy rozważyc w spos ˙ob następujący:

Z

b

a

f (x)dx ≃

n

X

i

=0

a

i

f (x

i

).

Rozważymy stosowane najczęściej w obliczeniach: metode trapezów i Simp-
sona.

17

18

ROZDZIAŁ 4. CAŁKOWANIE NUMERYCZNE

4.0.3

Metoda trapezów

Niech x

0

= a, x

1

= b, h = b − a ponownie wykorzystujemy twierdzenie z

rozdziału pierwszego

f (x) =

x − x

1

x

0

− x

1

f (x

0

) +

x − x

0

x

1

− x

0

f (x

1

)

Po podstawieniu:

Z

b

a

f (x)dx =

Z

x

1

x

0

x − x

1

x

0

− x

1

f (x

0

) +

x − x

0

x

1

− x

0

f (x

1

)

dx

+

1
2

Z

x

1

x

0

f ”(ζ(x))(x − x

0

)(x − x

1

)dx.

Wykorzystując twierdzenie o wartości średniej:

Z

x

1

x

0

f ”(ζ(x))(x − x

0

)(x − x

1

)dx = f ”(ζ)

Z

x

1

x

0

(x − x

0

)(x − x

1

)dx

= −

1
6

h

3

.

Ostatecznie po obliczeniach:

Z

b

a

f (x)dx =

h

2

(f (x

0

) + f (x

1

)) −

1

12

h

3

f ”(ζ).

Program w c/c++ oparty na tej metodzie możemy zrealizować w nastę-

pujący sposób, wynik uzyskamy tak jak w obliczeniach, oczywiście 0.5.

19

#include<iostream>
#include<cmath>
using namespace std;

double fun(double x)
{
return x;
}
int main()
{
double a=0.0;
double b=1.0;
int N=100;
double h=(b-a)/100;
double calka=0.;
double x0=a;
for(int i=0;i<N;i++)
{
calka=calka+(h/2.)*(fun(x0+i*h)+fun(x0+(i+1)*h));
}
cout<<calka<<endl;
return 0;
}

20

ROZDZIAŁ 4. CAŁKOWANIE NUMERYCZNE

4.0.4

Metoda Simpsona

Tym razem wybieramy trzy punkty : x

0

= a, x

2

= b, x

1

= a + h, h =

(b − a)/2

Z

b

a

f (x)dx =

Z

x

2

x

0

(x − x

1

)(x − x

2

)

(x

0

− x

1

)(x

0

− x

2

)

f (x

0

)dx

+

Z

x

2

x

0

(x − x

0

)(x − x

2

)

(x

1

− x

0

)(x

1

− x

2

)

f (x

1

)dx

+

Z

x

2

x

0

(x − x

0

)(x − x

1

)

(x

2

− x

0

)(x

2

− x

1

)

f (x

2

)dx

+

1
6

Z

x

2

x

0

f

(3)

(ζ(x))(x − x

0

)(x − x

1

)(x − x

2

)dx.

Po odpowiednim wycałkowaniu, uzyskujemy:

Z

x

2

x

0

f (x)dx =

h

3

(f (x

0

) + 4f (x

1

) + f (x

2

)) −

h

5

90

f

(4)

(ζ)

Kod programu w c/c++ realizujący powyższy problem:

#include<iostream>
#include<cmath>
using namespace std;

double fun(double x)
{

21

return x;
}
int main()
{
double a=0.0;
double b=1.0;
int N=100;
double h=(b-a)/100;
double calka=0.;
double x0=a;
for(int i=0;i<N/2;i++)
{
calka=calka+(h/3.)*(fun(x0+2*i*h)+4.*fun(x0+(2*i+1)*h)+fun(x0+(2*i+2)*h));
}
cout<<calka<<endl;
return 0;
}

22

ROZDZIAŁ 4. CAŁKOWANIE NUMERYCZNE

Rozdział 5

Numeryczne rozwiązywanie
równań

Kiedy mamy dane równanie kwadratowe typu

f (x) = ax

2

+ bx + c = 0,

bez problemu, przypominamy sobie odpowiednie wzory, które pozwalają nam
analitycznie wyznaczyć miejsca zerowe. Pojawia sie problem gdy mamy do
rozwiązania równanie postaci

f (x) = exp(−x) ln(x) − x ∗ x = 0.

Niestety już nie jesteśmy analitycznie rozwiązać to równanie. Możemy jedy-
nie podać przybliżona wartość miejsca zerowego.
Numeryczne metody rozwiązywanie równan f (x) = 0 z jedną niewiadoma:
1. Metoda połowienia (bisekcji)
Jeśli f ∈ C[a, b] i f (a) · f (b) < 0 to pierwiatek równania jest x

0

∈ [a, b].

Metoda oparta jest na procedurze zmiejszania przedziału w taki sposób że
po każdym kroku jest o połowę mniejszy.
przyklad
Rozważmy funkcje f (x) = x

3

+ 4x

2

− 10. Zadamy sobie pytanie czy w prze-

dziale < 1, 2 > znajduję się miejsce zerowe. Wyznaczamy:

f (1) = −5; f (2) = 6, f (1) · f (2) < 0.

Miejsce zerowe znajduję się w tym przedziale.

23

24

ROZDZIAŁ 5. NUMERYCZNE ROZWIĄZYWANIE RÓWNAŃ

Przykład
Numerycznie rozwiązujemy równanie typu f (x) = exp(−x) ln(x) − x ∗ x = 0..
Program w c/c++ realizujący ten problem może mieć postać:

#include<iostream>
#include<cmath>

using namespace std;

double fun(double x)
{
return exp(x)*log(x)-x*x;
}
int main()
{
double a=1.0;
double b=2.0;
double ab=b-a;
double delta=0.001;
double x0;
while( fabs(ab) > delta)
{
x0=(a+b)/2.0;
if(fun(a)*fun(x0) < 0.)
{

25

b=x0;
ab=b-a;
}
else
{
a=x0;
ab=b-a;
}
cout<<x0<<endl;
}
return 0;
}

i uzyskujemy rozwiązania kolejno krok po kroku:

1.5; 1.75; 1.625; 1.6875; 1.71875; 1.70312; 1.69531; 1.69141; 16336; 1.69434

2. Metoda Newtona

Jeśli f ∈ C

2

[a, b] tworzę ciąg typu

x

n

= x

n−

1

−

f (x

n−

1

)

f

′

(x

n−

1

)

.

Procedura trwa do momentu aż dla zadanej dokładności δ

|x

n

− x

n−

1

| ¬ δ.

26

ROZDZIAŁ 5. NUMERYCZNE ROZWIĄZYWANIE RÓWNAŃ

#include<iostream>
#include<cmath>

using namespace std;

double fun(double x)
{
return exp(x)*log(x)-x*x;
}

double fun1(double x)
{
return exp(x)*log(x)+exp(x)*(1./x)-2.*x;
}

int main()
{
double a=1.0;
double b=2.0;
double ab=b-a;
double delta=0.001;
double x0=(a+b)/2.;
double x1;
while( fabs(ab) > delta)
{
x1=x0-fun(x0)/fun1(x0);
ab=x1-x0;
x0=x1;
cout<<x1<<endl;
}
return 0;
}

Otrzymujemy w kolejnych krokach miejsc zerowe:1.7398; 1.69657; 1.6946; 1.6946.
3. Metoda siecznych

x

n

+1

= x

n

−

(x

n

− x

n−

1

)f

n

f

n

− f

n−

1

.

Rozdział 6

Równania rożniczkowe

Opis wielu zagadnień dla układów dynamicznych sprowadza się do równań
rożniczkowych pierwszego rzędu. Ogólnie możemy zapisać to w postaci:







































dy

1

dt

= g

1

(y1, y2, ..., y

n

, t)

dy

2

dt

= g

2

(y1, y2, ..., y

n

, t)

....

... ..............................

dy

n

dt

= g

n

(y1, y2, ..., y

n

, t)

Przyklad

Rozważmy oscylator harmoniczny w jednym wymiarze















dx

dt

= g

1

(x, v, t)

dv

dt

= g

2

(x, v, t)

gdzie g

1

(x, v, t) = v, g

2

(x, v, t) = −kx/m.

1. Metoda Eulera
Najprostszą metodą numerycznego przedstawienia równań różniczkowych jest
metoda Eulera. Polega ona na zastąpieniu pochodnych, pochodnymi nume-
rycznymi.

dy

dt

=

y

n

+1

− y

n

t

n

+1

− tn

=

y

n

+1

− y

n

h

= g(y

n

, t

n

)

co pozwala zapisać

y

n

+1

= y

n

+ hg

n

+ O(h

2

).

Algorytm oparty na metodzie Eulera jest mało dokładny i szybko robiega-
jący się. Dokładność metody można zwiększyć jedynie poprzez zwiększenie

27

28

ROZDZIAŁ 6. RÓWNANIA ROŻNICZKOWE

punktów początkowych. Dla wielu zagadnień nie jest to korzystne dlatego ze
zwykle zadany jest tylko jeden punkt zwany warunkiem początkowym.

2. Metoda Rungego-Kutty

Metoda Rungego-Kutty należy do tych algorytmow które nie wymagaj¸a zna-
jomości wielu punktów początkowych. Korzysta ona z dwóch rożnych szergów
Taylora dla zmiennych dynamicznych.

y(t + h) = y + hy

′

+

h

2

2

y“ +

h

3

3!

y

(3)

+ ...

= y + hg +

h

2

2

(g

t

+ gg

y

) +

h

3

3!

(g

tt

+ 2gg

tuy

+ g

2

g

yy

+ gg

2

y

+ g

t

g

y

) + ...

gdzie g

α,β

=

∂

2

g

∂α∂β

.

Równanie to możemy zapisać formalnie w postaci

y(t + h) = y(t) + α

1

k

1

+ α

2

k

2

+ ... + α

n

k

n

,

gdzie

k

n

= hg(y +

n−

1

X

l

=1

m

nl

k

l

+ h

n−

1

X

l

=1

m

nl

).

2a. Metoda Rungego-Kutty 2-go rzędu
W tym przypadku n=2, równania mają postać, wyższe rzędy zaniedbujemy,
traktując je jako zaniedbywalnie małe.

y(t + h) = y(t) + hg +

h

2

2

(g

t

+ gg

y

)











y(t + h) = y(t) + α

1

k

1

+ α

2

k

2

k

1

= hg(y, t),

k

2

= h(y + m

21

k

1

, t + m

21

h) = hg + m

21

h

2

(g

t

+ gg

y

)

wraz z ich postaciami po rozwinięciu w szerg Taylora z dokładnością do
wyrazów stopnia co najwyzej 2. Porównując oba równania uzyskujemy

y(t+h) = y(t)+hg +

h

2

2

(g

t

+gg

y

) ≡ y(t)+α

1

hg(y, t)+α

2

hg +m

21

h

2

(g

t

+gg

y

).

Poszukujemy współczynników spełniających relacje

(

1 = α

1

+ α

2

1
2

= α

2

m

21

Rozwiązania mogą mieć postać:

α

1

= α

2

=

1
2

, m

21

=

1
3

α

1

=

1
3

, α

2

=

2
3

, m

21

=

3
4

29

2b. Metoda Rungego-Kutty 4-go rzędu

Postępując podobnie jak w poprzedniej metodzie (w tym przypadku n=4)

y(t + h) = y(t) + α

1

k

1

+ α

2

k

2

+ α

3

k

3

+ α

4

k

4

,

gdzie i ich postać po rozwinięciu w szereg Taylora z dokładnością do wyrazów
stopnia co najwyżej 2

k

1

= hg(y, t),

k

2

= h(y + m

21

k

1

, t + m

21

h) = hg + m

21

h

2

(g

t

+ gg

y

),

k

3

= hg(y + m

31

k1 + m

32

k2, t + m

31

h + m

32

h),

k

4

= hg(y + m

41

h1 + m

32

k2 + m

33

k3, t + m

31

h + m

32

h + m

33

h).

Uzyskujemy następujące zależności pomiędzy współczynnikami

y(t + h) =

1
6

(k

1

+ 2k

2

+ 2k

3

+ k4),



















k

1

= hg(y, t)

k

2

= hg(y + 0.5k

1

, t + 0.5h)

k

3

= hg(y + 0.5k

3

, t + 0.5h)

k

4

= hg(y + k

3

, t + h)

Przyklad
Rozważmy oscylator harmoniczny z tłumieniem (y

1

= x, y

2

= v) z warunkami

początkowymi x(0) = 0, v(0) = v

0

. Równania mają postać:















dx

dt

= g

1

(x, v, t) = v

dv

dt

= g

2

(x, v, t) = −

kx

m

− bv

Wykorzystując metodę Rungego-Kutty 4-rzędu uzyskujemy

#include<iostream>
#include<cmath>

using namespace std;

double g1(double x, double v,double t)
{
return v;
}

30

ROZDZIAŁ 6. RÓWNANIA ROŻNICZKOWE

double g2(double x, double v, double t)
{
double a1=0.1;
double a2=0.1;
return -a1*x-a2*v;
}

int main()
{

double h=0.001;
int N=50;
double x0=1.;
double v0=0.;
double k11,k12,k13,k14;
double k21,k22,k23,k24;
double x,v;
double t;
t=0.;
x=x0;
v=v0;
for(int i=0;i<N;i++)
{
t=i*h;
k11=h*g1(x,v,t);
k21=h*g2(x,v,t);
k12=h*g1(x+k11/2.,v+k21/2.,t+h/2.);
k22=h*g2(x+k11/2.,v+k21/2.,t+h/2.);
k13=h*g1(x+k12/2.,v+k22/2.,t+h/2.);
k23=h*g2(x+k12/2.,v+k22/2.,t+h/2.);
k14=h*g1(x+k13,v+k23,t+h);
k24=h*g2(x+k13,v+k23,t+h);
x=x+(k11+2.*k12+2.*k13+k14)/6.;
v=v+(k21+2.*k22+2.*k23+k24)/6.;
cout<<t<<"\t"<<x<<endl;
cout<<t<<"\t"<<v<<endl;
}
return 0;
}

Bibliografia

[1] Tao Pang, Metody obliczeniowe w fizyce, PWN 2001

[2] Z.Fortuna, B.Macukow, J.Wąsowski, Metody Numeryczne, WNT 1998

[3] J.N.Bronsztajn, K.A. Siemiendiajew, Matematyka. Poradnik encyklope-

dyczny

, PWN 1999

[4] Jerzy Grębosz, Symfonia C++, Oficyna Kallimach, Kraków 1996

31

Skorowidz

aproksymacja

a. jednostajna 11
a. sredniokwadratowa 9
a. Pade’go 12
a. szeregami funkcyjnymi 11

interpolacja 3
całkowanie numeryczne

c. numeryczne trapezami 18
c. numeryczne Simpsona 20

numeryczne rozwiązywanie równań 23

metoda Newtona 23
metoda połowienia 25
metood siecznych 26

pochodna numeryczna

pierwsza pochodna numeryczna 2-

punktowa 15

pierwsza pochodna numeryczna 3-

punktowa 16

pierwsza pochodna numeryczna 5-

punktowa 16

równania różniczkowe 27

metoda Eulera 27
metoda Rungego-Kutty 28

32