METODY NUMERYCZNE

METODY NUMERYCZNE

Wykład 4.

Całkowanie numeryczne

Plan

• Wzór trapezów

• Złożony wzór trapezów

• Metoda ekstrapolacji Richardsona

• Metoda Romberga

• Metoda Simpsona – wzór parabol

• Metoda Gaussa

Całkowanie numeryczne - idea

b

a

dx ) x ( f I

X Y

f(x)

a b

b

a

dx x f( )

Całkę można przybliżyć sumą

n

i

f c

x

S

1

) (

x c

x

Kwadratury Newtona-Cotesa

Kwadratura Newtona – Cotesa należy do metod z

ustalonymi węzłami, polega na tym, że funkcja f(x) jest interpolowana wielomianem (np. Lagrange’a)

) x ( f )

x (

f

gdzie: _n

n n

n

n

( x ) a a x ... a x a x

f

Wówczas całka z funkcji f(x) może być przybliżana całką z funkcji interpolującej f

(x) n-tego stopnia

b a

n b

a

x f x

f

I ( ) ( )

Wzór trapezów

Zakładamy n = 1 czyli

f

( x ) a

a

x

dx x

a a

x f x

f I

b a b

a b

a

) (

) ( )

(

Szukamy współczynników a₀ i a₁

) 2 (

2 2

1 0

a a b

a b

a

Zakładamy, że prosta, która przybliża funkcję f(x) przechodzi przez punkty (a,f(a)) i (b,f(b)). Czyli:

a a a

a f a

f ( )

( )

a b

a b f b a

a f ( ) ( )

0

Wzór trapezów

trapezu pole

dx x

f

b

a

) (

X Y

f(x)

a b

f_I(x)

2

) ( )

) ( (

)

( f a f b

a b

dx x

f

b a

f(a) f(b)

Wzór trapezów

t t t

v 9 . 8

2100 140000

140000 ln

2000 )

(

Przykład 1:

Prędkość rakiety w przedziale czasu od t₁=8 s do t₂=30 s jest opisana wzorem:

a) Przy pomocy metody trapezów znajdź przemieszczenie rakiety w tym przedziale czasu czyli Δx=x(t₂)-x(t₁)

b) Wyznacz błąd względny odniesiony do wartości dokładnej (ang.

true relative error)

Wzór trapezów

2

) ( )

) (

( f a f b

a b

I a 8 s b 30 s

t t t

f 9 . 8

2100 140000

140000 ln

2000 )

(

) 8 ( 8 . ) 9

8 ( 2100 140000

140000 ln

2000 )

8 ( f

) 30 ( 8 . ) 9

30 ( 2100 140000

140000 ln

2000 )

30 ( f

s m / 27

. 177

s m / 67

. 901 a)

2

67 . 901 27

. ) 177 8

30 (

I 11868 m

Wzór trapezów

30 8

8 . 2100 9

140000

140000 ln

2000 t dt

x t 11061 m

b)

11061 100

11868 11061

t

7 .2959 %

Błąd względny:

Złożony wzór trapezów

n a h b

Błąd, jak pokazuje poprzedni przykład, jest zbyt duży. Można zaproponować podział przedziału całkowania na n segmentów,

każdy o długości h.

h a

h a

h a

h a

h a

h a h

a a b

a

dx x f dx

x f dx

x f dx

x f dx

x f I

2 3

2

4

3

) ( )

( )

( )

( )

(

dla n=4 ^X

Y

f(x)

a b

4 a a b

2b4a

a 3b4a

a

b

dx ) x (

f _{f ( a ) f ( a ih ) f ( b )}

n a

b

2 2

b

h ) n ( a h

) n ( a

h ) n ( a h

a

h a h

a

a

dx ) x ( f dx

) x ( f ...

dx ) x ( f dx

) x ( f

1 1

2 b 2

a

dx ) x ( f

2 ...

) 2 (

) (

2

) (

)

( f a h f a h

h h a

f a

h f

2

) ( )

) 1 (

] ( ) 1 (

( [

... f a n h f b

h n

a b

Złożony wzór trapezów

t t t

v 9 . 8

2100 140000

140000 ln

2000 )

(

Przykład 2:

Prędkość rakiety w przedziale czasu od t₁=8 s do t₂=30 s jest opisana wzorem:

a) Przy pomocy metody trapezów znajdź przemieszczenie rakiety w tym przedziale czasu czyli Δx=x(t₂)-x(t₁) przyjmując n=2

b) Wyznacz błąd względny odniesiony do wartości dokładnej (true relative error)

Złożony wzór trapezów

a) ( ) 2 ( ) ( ) 2

1 1

b f ih

a f a

n f a

I b

i

n = 2 a = 8 s b = 30 s

s

n a

h b 11

2 8 30

) 30 ( )

( 2

) 8 ) (

2 ( 2

8

30

1

f ih

a f f

I

i

) 30 ( )

19 ( 2 ) 8 4 (

22 f f f 177 27 2 484 75 901 67

4

22 . ( . ) .

Złożony wzór trapezów

b)

30

8

8 2100 9

140000

140000

2000 . t dt

ln t

x 11061 m

wartość dokładna wynosi:

Błąd względny to

:

11061 100

11266 11061

t

1 . 8534 %

Złożony wzór trapezów

n Δx E_t

1 11868 -807 7.296 ---

2 11266 -205 1.853 5.343

3 11153 -91.4 0.8265 1.019 4 11113 -51.5 0.4655 0.3594 5 11094 -33.0 0.2981 0.1669 6 11084 -22.9 0.2070 0.09082 7 11078 -16.8 0.1521 0.05482 8 11074 -12.9 0.1165 0.03560

t % _a %

Złożony wzór trapezów

Błąd bezwzględny metody z pojedynczym segmentem

b a

), (

"

) f a b

E

(

12

3

gdzie jest punktem w a,b

Błąd w metodzie złożonej (wielosegmentowej) jest sumą błędów dla każdego segmentu. Błąd pojedynczego segmentu wynosi:

h a

a ), (

"

a f )

h a

E (

3

1

12

) (

"

h f

1 3

12

Analiza błędu dla wzoru trapezów

Podobnie:

ih a

h i

a h f

i a ih

E_i a "( _i), ( 1) _i

12

) ) 1 ( (

)

( ³

) (

"

h f 12

3

dla n:

b h

) n

( a ), (

"

h f ) n

( a

E

b

1

12

1

) (

"

h f

n 3

Analiza błędu dla wzoru trapezów

Całkowity błąd w metodzie złożonej jest sumą błędów pojedynczych segmentów:

n i

i

t

E

E

1

n i

i

) (

"

h f

1 3

12 n

) (

"

f n

) a b

(

n i

i 1

2 3

12

Wyrażenie

n

) (

"

f

n i

i 1

jest przybliżoną średnią wartością drugiej pochodnej w przedziale

a x b

2

1 E

n

Analiza błędu dla wzoru trapezów

Poniższa tablica dla całki ³⁰

8

8 2100 9

140000 140000

2000 . t dt

ln t

w funkcji liczby segmentów n. Widać, że gdy liczba segmentów jest podwajana, błąd E_t zmniejsza się w przybliżeniu czterokrotnie.

Et _t % _a %

n Value

2 11266 -205 1.854 5.343

4 11113 -51.5 0.4655 0.3594

8 11074 -12.9 0.1165 0.03560

Analiza błędu dla wzoru trapezów

Całkowanie metodą Romberga

Metoda Romberga jest rozszerzeniem metody

trapezów i daje lepsze przybliżenie całki poprzez

zasadniczą redukcję błędu (true error)

n

E

C

n

t

TV I

E

Ekstrapolacja Richardsona

Błąd (E_t true error) w n-segmentowym wzorze trapezów wynosi

gdzie C jest współczynnikiem proporcjonalności Stąd:

wartość dokładna (true

wartość przybliżona np.

I

TV n

C

2 2

2

gdy segment zostaje zmniejszony o połowę

3

2 2

n n

n

I I I

TV

Ekstrapolacja Richardsona

Można pokazać, że

Ze wzorów: TV I

n C

2

I

TV n

C

2 2

2

otrzymujemy:

t t t

v 9 . 8

2100 140000

140000 ln

2000 )

(

Przykład 3:

Prędkość rakiety w przedziale czasu od t₁=8 s do t₂=30 s jest opisana wzorem:

a) Przy pomocy ekstrapolacji Richardsona znajdź przemieszczenie rakiety w tym przedziale czasu czyli Δx=x(t₂)-x(t₁)

b) Wyznacz błąd względny odniesiony do wartości dokładnej (true relative error)

Ekstrapolacja Richardsona

t % _a %

n Δx E_t

1 11868 -807 7.296 ---

2 11266 -205 1.853 5.343

3 11153 -91.4 0.8265 1.019 4 11113 -51.5 0.4655 0.3594 5 11094 -33.0 0.2981 0.1669 6 11084 -22.9 0.2070 0.09082 7 11078 -16.8 0.1521 0.05482 8 11074 -12.9 0.1165 0.03560

Tabela wyników ze złożonego wzoru trapezów do n=8 segmentów

a)

m

I

11266 I

11113 m

3

2 2

n n

n

I I I

TV dla n=2

3

2 4

4

I I I

TV 3

11266 11113

11113 m 11062

Ekstrapolacja Richardsona

b)

30 8

8 2100 9

140000

140000

2000 . t dt

ln t

x 11061 m

Stąd

11062 11061

E

^{1 m}

Ekstrapolacja Richardsona

Wartość dokładna to:

.

c) Błąd względny

11061 100

11062 11061

t

0 .00904 %

Ekstrapolacja Richardsona

n Δx (m)

wzór

trapezów wzór trapezów

Δx (m)

ekstrapolacja

Richardsona ekstrapolacja Richardsona

12 4

11868 11266 11113

7.296 1.854 0.4655

11065-- 11062

0.03616-- 0.009041

Porównanie wyników z metodą złożoną trapezów

t % _t %

Całkowanie metodą Romberga stosuje ten sam wzór co ekstrapolacja Richardsona. Jednakże, metoda

Romberga jest to algorytm rekurencyjny.

Przypomnijmy:

3

2 2

n n

n

I I I

TV Można zapisać

3

2 2

2

n n

R n n

I I I

I 4

1

2 2

n n

n

I I I

Metoda Romberga

Wartość dokładna TV jest zastąpiona przez wynik całkowania metodą Richardsona

Znak jest zastąpiony przez znak równości.

n R

I

Estymowana wartość dokładna wynosi:

2

Ch

I

TV

gdzie Ch

jest wartością błędu przybliżenia Metoda Romberga

Następna wartość całki (podwajając liczbę segmentów) wynosi:

3

2 4

4 4

n n

R n n

I I I

I

Estymowana wartość dokładna wynosi teraz:

15

2 4

4 n R n R

n R

I I I

TV

1 2 4

1 1

1 1

1

I I , k

I

I

Wskaźnik k reprezentuje rząd ekstrapolacji

k=1 odpowiada wartościom uzyskanym ze wzoru trapezów Ogólne wyrażenie w metodzie Romberga

Metoda Romberga

k=2 odpowiada wartościom uzyskanym z błędem O(h

)

Wskaźnik j reprezentuje dokładność; j+1 daje całkę

wyznaczoną dokładniej niż j

t t t

v 9 . 8

2100 140000

140000 ln

2000 )

(

Przykład 4:

Prędkość rakiety w przedziale czasu od t₁=8 s do t₂=30 s jest opisana wzorem:

a) Przy pomocy wzoru Romberga znajdź przemieszczenie rakiety w tym przedziale czasu czyli Δx=x(t₂)-x(t₁)

b) Wyznacz błąd względny odniesiony do wartości dokładnej (true relative error)

Metoda Romberga

t % _a %

Tabela wyników ze złożonego wzoru trapezów do n=8 segmentów

n Δx E_t

1 11868 -807 7.296 ---

2 11266 -205 1.853 5.343

3 11153 -91.4 0.8265 1.019 4 11113 -51.5 0.4655 0.3594 5 11094 -33.0 0.2981 0.1669 6 11084 -22.9 0.2070 0.09082 7 11078 -16.8 0.1521 0.05482 8 11074 -12.9 0.1165 0.03560

Na podstawie tabeli odczytujemy:

11868

1

I

I

11266

11113

3

I

I

11074

Metoda Romberga

W pierwszym przybliżeniu:

3

1 1 2

1 2

1 1

2

, ,

, ,

I I I

I

3

11868 11266

11266

Podobnie,

3

2 , 1 3

, 1 3

, 1 2

, 2

I I I

I

3

11266 11113

11113

11062

3

3 , 1 4

, 1 4

, 1 3

, 2

I I I

I

3

11113 11074

11074

11061

Metoda Romberga

W drugim przybliżeniu,

15

1 2 2

2 2

2 1

3

, ,

, ,

I I I

I

15

11065 11062

11062 11062 Podobnie,

15

2 , 2 3

, 2 3

, 2 2

, 3

I I I

I

11062 11061

Metoda Romberga

Dla trzeciego rzędu,

63

1 3 2

3 2

3 1

4

, ,

, ,

I I I

I

63

11062 11061

11061

m 11061

Metoda Romberga

11868

11262

11113

11074

11065

11062

11061

11062

11061

11061 1-segment

2-segment

4-segment

8-segment

Rząd 1 Rząd 2 Rząd 3

Poprawione wartości całki metodą Romberga

Metoda Romberga

Metoda Simpsona

Metoda trapezów była oparta na przybliżeniu funkcji podcałkowej f(x) wielomianem stopnia pierwszego.

W metodzie Simpsona wykorzystuje się wielomiany stopnia drugiego, jest to tzw. metoda parabol.

b

a b

a

dx ) x ( f dx

) x ( f

I

gdzie:

2 2 1

0

2

( x ) a a x a x

f

X Y

f(x)

a b

f_I(x)

Metoda Simpsona

Równanie paraboli dla 3 punktów:

)), a ( f , a (

b , f a

b , a

2 2

)) b ( f , b (

2 2 1

0

2

( a ) a a a a a

f )

a ( f

2 2

1 0

2

2 2 2

2

b a a

b a a

b a f a

b

f a

Metoda Simpsona

2 2

2 2

0

2

4 2

b ab

a

) a ( f b ) a ( b abf

abf a )

b ( abf )

b ( f a a

2 1 2

2

4 2 3

2 3 4

b ab

a

) b ( b bf

bf a )

a ( bf )

b ( b af

af a )

a ( af a

2 2 2

2 2 2 2

b ab

a

) b ( b f

f a )

a ( f a

Wyznaczone współczynniki funkcji f

(x) to:

Metoda Simpsona

b

a

dx ) x ( f

I

b

a

dx x

a x

a

a

b

a

a x a x

x

a 2 3

3 2 2

1 0

3 2

3 3

2 2

2 1 0

a a b

a a b

) a b

(

a

Wówczas:

Metoda Simpsona

) b ( b f

f a )

a ( a f

dx b ) x ( f

b

a

4 2

2

6

2 a h b

Co daje:

) 2 (

4 )

3 ( )

2

( a b f b

f a

h f dx

x f

b

a

wzór parabol

Metoda Simpsona w wersji złożonej

) ...

x ( f ) x ( f )

x ( ) f

x x

( dx ) x ( f

b

a

6

4

0 0

2

) ...

x ( f ) x ( f )

x ( ) f

x x

( 6

4

2 2

4

f(x)

. . .

x₀ x₂ x_n-2 x_n x

...

dx ) x ( f dx

) x ( f dx

) x ( f

x

x x

x b

a

4

2 2

0

n

n n

n

x

x x

x

dx ) x ( f dx

) x ( f ....

2 2

4

) ...

x ( f ) x

( f )

x ( ) f

x x

(

...

6

4

4 4

2

h x

x

2

n

...,

,

,

i 2 4

Metoda Simpsona

) ...

x ( f ) x ( f )

x ( h f

dx ) x ( f

b

a

6

2

4

) ...

x ( f ) x ( f )

x ( h f

6

2

4

) ...

x ( f ) x

( f )

x (

h f

6

2

4

6

2 f ( x

) 4 f ( x

) f ( x )

h

Metoda Simpsona

b

a

dx ) x (

f ^{h f ( x}

^{) 4 f ( x}

^{) f ( x}

^{) ... f ( x}

^{) ...}

3

)}]

( )

( ...

) ( )

( 2

... f x

f x

f x

f x

) ( )

( 2

) ( 4

) 3 (

2 2 1

1

0 n

n

even i

i

i n

odd i

i

i

f x f x

x f x

h f

) ( )

( 2

) ( 4

) 3 (

2 1

0 n

n

i n

i

f x f x

x f x

n f

a

b

Wartości przybliżone

przykładu stosując regułę 1/3 Simpsona z wieloma segmentami

n Wartość przybliżona E_t ^{|Єt |} 24

68 10

11065.72 11061.64 11061.40 11061.35 11061.34

4.380.30 0.060.01 0.00

0.0396%

0.0027%

0.0005%

0.0001%

0.0000%

Metoda Simpsona – analiza błędu

Błąd dla jednego segmentu

b a f a b

E

( ),

2880 )

(

) (

) f x x

E (

5 0 2

1

2880

), (

h f

1 4

5

90 )

( ) f

x x

E (

5 2 4

2

2880 h f

( ),

2 4

5

90

2 1

0

x

x

4 2

2

x

x

Metoda Simpsona – analiza błędu

Błąd w metodzie wielosegmentowej

) x

x

(

Całkowity błąd

2 1 n

i

i

t

E

E

2 1

) 4 ( 5

) 90 (

n

i

f

h

1

) 4 ( 5

5

) 90 (

) (

n

i

f

n a b

n f n

a b

n

i i

2 1

) 4 (

4

5

( )

90 ) (

Metoda Simpsona – analiza błędu

) 4 ( 4

5

90 )

( f

n a E

b

średnia wartość pochodnej

n f f

n

i

i 2

1

) 4 ( )

4

(

( )

Metoda Gaussa

b

a

dx ) x ( f

I c

f ( x

) c

f ( x

)

Całkę metodą kwadratury Gaussa przedstawia wzór:

Punkty x₁ i x₂ , w których określamy wartość funkcji podcałkowej nie są ustalone (jak poprzednio na granicach przedziału czyli a i b), ale

stałe współczynniki

Metoda Gaussa

. x a x

a x

a a

) x (

f

b

a b

a

dx x

a x

a x

a a

dx ) x (

f

b

a

a x a x

a x x

a 2 3 4

4 3 3

2 2

1 0

4 3

2

4 4

3 3

3 2

2 2

1 0

a a b

a a b

a a b

a b

a

Niewiadome x

, x

, c

, c

znajduje się zakładając, że

funkcja podcałkowa spełnia warunek:

Metoda Gaussa

3 2 3 2

2 2 2

1 0

2 3

1 3 2

1 2 1

1 0

1

a a x a x a x c a a x a x a x

c dx

) x ( f

b

a

3 2 2 3

1 1 3 2

2 2 2

1 1 2 2

2 1

1 1 2

1

0

c c a c x c x a c x c x a c x c x

a

b

a

dx ) x ( f

I c

f ( x

) c

f ( x

)

b

a b

a

dx x

a x

a x

a a

dx ) x (

f

4 4

3 3

2

2

a b a b a

b

z jednej strony

z drugiej strony

Metoda Gaussa

2

1

c

c a

b

2 2 1

1 2

2

2 a c x c x b

2 2 2 2

1 1 3

3

3 a c x c x

b

2 2 3

1 1 4

4

4 a c x c x b

3 2 2 3

1 1 3 2

2 2 2

1 1 2 2

2 1

1 1 2

1

0

c c a c x c x a c x c x a c x c x

a

4 3

2

4 4

3 3

3 2

2 2

1 0

a a b

a a b

a a b

a b

a

Metoda Gaussa

1

1

a b

a x b

2 3

1

2

2

a b

a x b

Rozwiązując układ równań

:

2 1

c c

a b

2 2 1

1 2

2

2 a c x c x b

2 2 2 2

1 1 3

3

3 a c x c x b

3 2 2 3

1 1 4

4

4 a c x c x b

otrzymujemy dla metody dwupunktowej

:

Metoda Gaussa

3 2 1 2

2 3 2

1 2

2 )

(

a b

a f b

a b

a b

a f b

a b

x f c x

f c dx

x f

b

a

W dwu-punktowej metodzie Gaussa, całka funkcji f(x) wyraża się wzorem:

) x ( f c .

. . . . . . )

x ( f c )

x ( f c dx

) x (

f

b

a

2 2

1 1

Uogólniając dla n punktów:

1

1 1

n

i

c

g ( x

) dx

) x ( g

n współczynniki argumenty funkcji 2 c₁ = 1.000000000

c₂ = 1.000000000 x₁ = -0.577350269 x₂ = 0.577350269 3 c₁ = 0.555555556

c₂ = 0.888888889 c₃ = 0.555555556

x₁ = -0.774596669 x₂ = 0.000000000 x₃ = 0.774596669 4 c₁ = 0.347854845

c₂ = 0.652145155 x₁ = -0.861136312 x₂ = -0.339981044

Metoda Gaussa

W n-punktowej metodzie Gaussa, współczynniki c_i oraz argumenty x_i są stabelaryzowane dla całek w granicach od -1 do 1

n współczynniki argumenty funkcji

5 c₁ = 0.236926885 c₂ = 0.478628670 c₃ = 0.568888889 c₄ = 0.478628670 c₅ = 0.236926885

x₁ = -0.906179846 x₂ = -0.538469310 x₃ = 0.000000000 x₄ = 0.538469310 x₅ = 0.906179846 6 c₁ = 0.171324492

c₂ = 0.360761573 c₃ = 0.467913935 c₄ = 0.467913935 c₅ = 0.360761573 c₆ = 0.171324492

x₁ = -0.932469514 x₂ = -0.661209386 x₃ = -0.2386191860 x₄ = 0.2386191860 x₅ = 0.661209386 x₆ = 0.932469514

Metoda Gaussa

Jeżeli dane są w tablicach dla

1 1

dx ) x (

g

b

a

dx ) x (

f

Granice całkowania

a , b

1, 1

Niech

x mt c

Dla

x a , _{t 1}

, b

Dla

x t 1

Wynika stąd, że:

2 a m b

Metoda Gaussa

Stąd:

2 2

a t b

a x b

a dt dx b

2

Ostatecznie:

1

1

2 2 2

)

( b a b a dt

a t f b

dx x

f

b

a

Metoda Gaussa

t t t

v 9 . 8

2100 140000

140000 ln

2000 )

(

Przykład 5:

Prędkość rakiety w przedziale czasu od t₁=8 s do t₂=30 s jest opisana wzorem:

a) Przy pomocy dwu-punktowej metody Gaussa znajdź

przemieszczenie rakiety w tym przedziale czasu czyli Δx=x(t₂)- x(t₁)

b) Wyznacz błąd względny odniesiony do wartości dokładnej (true

Metoda Gaussa

Najpierw zmieniamy granice całkowania z [8,30] na [-1,1]

1 1 30

8

2

8 30 2

8 30 2

8

30 f x dx

dt ) t ( f

1 1

19 11

11 f x dx

Metoda Gaussa

Następnie odczytujemy z tablic dla n=2 c₁ 1.000000000 577350269

1 0.

x

000000000

2 1.

c

577350269

2 0.

x

Korzystamy ze wzoru kwadratury Gaussa

19 11

11 19

11 11

19 11

11

1 1

x f

c x

f c dx

x f

19 5773503

0 11 11

19 5773503

0 11

11 f ( . ) f ( . )

) .

( f )

. (

f 12 64915 11 25 35085 11

) .

( )

.

( 296 8317 11 708 4811 11

m .44 11058

Metoda Gaussa

Skorzystaliśmy z tego, że:

) .

( ) .

. ln (

) .

(

f 9 8 12 64915

64915 12

2100 140000

140000 2000

64915 12

8317 296.

) .

( ) .

. ln (

) .

(

f 9 8 25 35085

35085 25

2100 140000

140000 2000

35085 25

4811 708.

Metoda Gaussa

Błąd względny, _t

. %

. .

t

100

34 11061

44 11058 34

11061

c)

Błąd bezwzględny (true error)

b)

E

44 . 11058 34

. 11061 E

m 9000 .

2

Metoda Gaussa