Mariusz PYRZ

SIMR (PW), Instytut Pojazdów

Metody numeryczne w mechanice

Całkowanie numeryczne

8.

Przybli

ż

one obliczanie całek oznaczonych

Aby obliczyć wartość całki oznaczonej

gdzie f: R->R spełnia niezbędne warunki różniczkowania i całkowania

przybliżamy funkcję f(x) w przedziale [a,b] za pomocą wielomianu P

N

(x)

łatwego do scałkowania i na jego podstawie obliczamy przybliżoną wartość I

a

( )

b

b

I

f x dx

=

∫

MNM Całkowanie numeryczne M.Pyrz 04-2014

2

Aby uzyskać dokładniejsze przybliżenie rozwiązania
(ale uniknąć zwiększania stopnia wielomianu interpolacyjnego P

N

(x) w [a,b]),

dzielimy przedział [a,b] na podprzedziały, w każdym z nich obliczamy

oddzielnie wartość całki i następnie sumujemy wszystkie oszacowania.

( )

( )

b

b

a

N

b

b

I

P x dx

f x dx

=

≈

∫

∫

Obliczanie całki w przedziale [a,b]

Zastosowanie:
obliczanie powierzchni, objętości obszarów o znanym brzegu, charakterystyk przekrojów

(środek ciężkości, moment bezwładności), itp.

3

Rys. 8.1

(

) /

h

b

a

N

= −

MNM Całkowanie numeryczne M.Pyrz 04-2014

Metoda prostok

ą

tów

Oszacowaniem wartości całki w przedziale [a,b] jest pole prostokąta

( )

( )(

)

b

a

a

f x dx

I

f a b

a

≈

=

−

∫

( )(

)

a

I

f b b

a

=

−

f(a) – lewa granica przedziału

f(b) - prawa granica przedziału

4

Rys. 8.2

(

)(

)

2

a

a

b

I

f

b

a

+

=

−

f((a+b)/2) – środek przedziału

MNM Całkowanie numeryczne M.Pyrz 04-2014

Numeryczne obliczanie całki - kwadratura

Dla skończonej liczby punktów (węzłów)

można wyznaczyć wielomian interpolacyjny (przechodzący przez te węzły)

Wartość całki (dla wielomianu) równa jest sumie ważonej wartości całkowanej

0

1

2

...

N

a

x

x

x

x

b

=

< <

< <

=

[ , ]

0,1,...,

k

x

a b

k

N

∈

=

0

( )

( ) ( )

N

N

i

i

P x

L x f x

=

∑

5

funkcji liczonej w kilku odpowiednich punktach (nazywanej także kwadraturą)

A

k

są stałymi niezależnymi od f(x):

( )

b

k

k

b

A

L x dx

=

∫

0

ˆ

( )

(

),

[ , ]

b

N

N

N

k

k

k

k

a

I

P x dx

A f x

x

a b

=

=

=

∈

∑

∫

MNM Całkowanie numeryczne M.Pyrz 04-2014

Dokładno

ść

w obliczaniu całek

Współczynniki A

k

i węzły kwadratury x

k

powinny być tak dobrane,

aby zminimalizować błąd aproksymacji

ˆ

N

N

R

I

I

=

−

Aby kwadratura była zbieżna, musi zachodzić

0

N

N

R

→∞

→

6

MNM Całkowanie numeryczne M.Pyrz 04-2014

Obliczanie przybli

ż

onej warto

ś

ci całek

Metody Newtona-Cotesa

wykorzystują interpolację funkcji f(x) za pomocą

wielomianów Lagranga L

N

, budowanych na węzłach równoodległych

Metoda trapezów (najprostsza z metod Newtona-Cotesa) dokonuje interpolacji funkcji f(x) za

0

0,

( )

(

)

N

N

j

N

n

n

j

j n

n

j

x

x

L

x

f x

x

x

=

=

≠

−

=

−

∑

∏

,

(

) /

0,1,...,

n

x

a

nh

h

b a

N

n

N

= +

= −

=

7

Metoda trapezów (najprostsza z metod Newtona-Cotesa) dokonuje interpolacji funkcji f(x) za
pomocą wielomianu pierwszego stopnia (tj. odcinka prostej).

Metody Gaussa

stosują interpolacje budowane na węzłach

ξ

i

rozłożonych

nierównomiernie i dobieranych tak, aby anulować błąd aproksymacji R

N

.

MNM Całkowanie numeryczne M.Pyrz 04-2014

1

1

1

( )

( )

n

i

i

i

I

f x dx

f

w

ξ

+

=

−

=

=

∑

∫

Metoda trapezów

Przedział całkowania [a,b] dzielony jest na N

podprzedziałów o długości h=(b-a)/N.

W każdym z podprzedziałów krzywa f(x) jest

zastępowana odcinkiem prostej.

Otrzymana aproksymacja jest rzędu O(h

2

) co

oznacza, że różnica między dokładną wartością całki
a wartością obliczoną jest proporcjonalna do h

2

).

8

a wartością obliczoną jest proporcjonalna do h

2

).

f x dx

x

x

i

i

( )

+

z

1

[

]

1

ˆ

(

)

( )

2

i

i

i

h

I

f x

f x

−

=

+

MNM Całkowanie numeryczne M.Pyrz 04-2014

Algorytm

Ustalić h=(b-a)/N oraz x

i

=a+ih (mamy w sumie N+1 punktów x

i

, x

0

=a, x

N

=b ).

W podprzedziale [x

i

, x

i+1

] zastępujemy wartość całki

polem trapezu

Metoda trapezów

1

1

1

ˆ

( )

( )

( )

2

(

)

2

b

N

N

i

N

N

i

i

a

h

f x dx

I

R

f a

f b

f a

ih

R

−

=

=





=

+

=

+

+

+

+









∑

∑

∫

1

2

1

ˆ

ˆ

ˆ

ˆ

( )

...

b

N

N

i

i

a

f x dx

I

I

I

I

=

≈ + + +

=

∑

∫

9

Błąd

3

[ , ]

max

( )

12

N

x

a b

h

R

N

f

x

∈

′′

≤

MNM Całkowanie numeryczne M.Pyrz 04-2014

Metoda Simpsona

Przedział całkowania [a,b] dzielony jest na N
podprzedziałów o długości h=(b-a)/N.
W każdym podprzedziale krzywa f(x) jest
zastępowana wielomianem interpolującym drugiego
stopnia (czyli łukiem paraboli).
Uzyskana aproksymacja jest rzędu O(h

4

) (czyli

różnica między wartością dokładną całki a wartością
obliczoną jest proporcjonalna do h

4

).

10

obliczoną jest proporcjonalna do h

4

).

2

2

2

( )

i

i

x

x

f x dx

+

∫

[

]

2

2

2

1

2

2

ˆ

(

)

4 (

)

(

)

3

i

i

i

i

h

I

f x

f x

f x

+

+

=

+

+

MNM Całkowanie numeryczne M.Pyrz 04-2014

Algorytm

Ustalić h=(b-a)/2N oraz x

i

=a+ih, (mamy teraz 2N+1 punktów x

i

, x

0

=a, x

2N

=b ).

W podprzedziale [x

2i

, x

2i+2

] zastępujemy wartość całki

polem „pod wielomianem”

Metoda Simpsona

1

/ 2 1

/ 2 1

2

0

1

0

ˆ

( )

( )

( )

2

(

2 )

4

(

(2

1) )

3

b

N

N

N

i

N

N

i

i

i

a

h

f x dx

I

R

f a

f b

f a

ih

f a

i

h

R

−

−

−

=

=

=





=

+

==

+

+

+

+

+

+

+









∑

∑

∑

∫

1

0

2

4

2(

1)

2

0

ˆ

ˆ

ˆ

ˆ

ˆ

( )

...

b

N

N

i

i

a

f x dx

I

I

I

I

I

−

−

=

≈ + + + +

=

∑

∫

11

5

(4)

[ , ]

max

( )

180

N

x

a b

h

R

N

f

x

∈

≤

MNM Całkowanie numeryczne M.Pyrz 04-2014

Błąd

Metoda Boole’a

Przedział [a,b] dzielony jest na podprzedziały. W każdym z nich

funkcja f(x) jest zastępowana wielomianem 4 stopnia.

Uzyskana aproksymacja jest rzędu O(h

6

)

Algorytm

Ustalić h=(b-a)/4N, x

i

=a+ih.

W każdym podprzedziale [x

4i

, x

4i+4

] wartość całki

f x dx

x

i

( )

4

4

+

z

12

W każdym podprzedziale [x

4i

, x

4i+4

] wartość całki

jest przybliżana wzorem

f x dx

x

i

( )

4

z

[

]

4

4

4

1

4

2

4

3

4

4

2

ˆ

7 (

) 32 (

) 12 (

) 32 (

) 7 (

)

45

i

i

i

i

i

i

h

I

f x

f x

f x

f x

f x

+

+

+

+

=

+

+

+

+

1

4

0

ˆ

( )

b

N

i

i

a

f x dx

I

−

=

≈

∑

∫

MNM Całkowanie numeryczne M.Pyrz 04-2014

Oszacowanie bł

ę

du

Oszacowanie wartości błędu całki w przedziale [a,b]:

Dla [a,b] podzielonego na N podprzedziałów:

(błąd liczymy rozwijając f(x) w szereg Taylora, całkujemy żeby obliczyć I, odejmujemy I

a

,

wynik wyrażamy w funkcji h – szerokości podprzedziału)

a

E

I

I

= −

1

N

i

i

E

E

=

=

∑

13

Metoda prostokątów

Metoda prostokątów (f w środku przedz.)

Metoda trapezów

Metoda Simpsona

2

( )

2

b

a

E

f h

O h

−

′

=

=

2

2

(

)

24

b

a

E

f h

O h

−

′′

=

=

2

2

(

)

12

b

a

E

f h

O h

−

′′

=

=

4

4

(

)

180

IV

b

a

E

f

h

O h

−

=

=

MNM Całkowanie numeryczne M.Pyrz 04-2014

Metoda Romberga

Metoda ta pozwala na ustalenie dokładności wyniku i prowadzenie obliczeń

iteracyjnych aż do osiągnięcia żądanej dokładności. Algorytm metody

bazuje na oryginalnym wykorzystaniu metody trapezów.

•T

0

(h)

– oszacowanie całki uzyskane na podstawie metody trapezów o kroku h

•T

0

(h/2)

– dokładniejsze oszacowanie wartości całki I otrzymane dla kroku h/2

•T (h/4)

– oszacowanie całki I uzyskane przy kroku h/4

14

•T

0

(h/4)

– oszacowanie całki I uzyskane przy kroku h/4

•T

0

(h/8) …

•T

0

(h/2

n

)

- ciąg jest zbieżny do dokładnej wartości I.

Można wykazać, że można utworzyć drugi ciąg T

1

(h/2

n-1

) w którym element

T

1

(h/2

k

) stanowi lepsze przybliżenie całki niż T

0

(h/2

k

).

MNM Całkowanie numeryczne M.Pyrz 04-2014

Metoda Romberga

Otrzymuje się

W analogiczny sposób można otrzymać trzeci ciąg

h

h

FG IJ FG IJ

T

h

T

h

T

h

k

n

k

k

k

1

0

1

0

2

4

2

2

4 1

0

1

F

HG

I

KJ

=

F

HG

I

KJ

−

F

HG

I

KJ

−

∈

−

+

[ ,

]

15

W sformułowaniu ogólnym

T

h

T

h

T

h

k

n

k

k

k

2

2

1

1

1

2

2

4

2

2

4

1

0

2

F

HG

I

KJ

=

F

HG

I

KJ

−

F

HG

I

KJ

−

∈

−

+

[ ,

]

T

h

T

h

T

h

k

n

p

p

k

p

p

k

p

k

p

2

4

2

2

4

1

0

1

1

1

F

HG

I

KJ

=

F

HG

I

KJ

−

F

HG

I

KJ

−

∈

−

−

+

−

[ ,

]

MNM Całkowanie numeryczne M.Pyrz 04-2014

Metoda Romberga

Grupujemy elementy w następującej tablicy

T h

T

h

T

h

T

h

T

h

T h

T

h

T

h

T

h

T h

T

h

T

h

n

n

n

n

n

n

0

0

0

2

0

1

0

1

1

1

2

1

1

2

2

2

2

2

2

2

2

2

2

2

2

2

b g

b g

b g

F

HG

I

KJ

F

HG

I

KJ

F

HG

I

KJ

F

HG

I

KJ

F

HG

I

KJ

F

HG

I

KJ

F

HG

I

KJ

F

HG

I

KJ

F

HG

I

KJ

−

−

−

−

−

⋯

⋯

⋯

16

Można wykazać, że ciąg T

k

(h) zbiega się szybciej do I niż ciąg T

0

(h/2

k

).

T

h

T

h

T h

n

n

n

n

2

2

2

2

1

1

2

2

2

b g

b g

b g

HG KJ

HG KJ

F

HG

I

KJ

−

−

−

⋯

⋯

MNM Całkowanie numeryczne M.Pyrz 04-2014

Metoda Romberga

Algorytm

Obliczyć metodą trapezów T

0

(h), T

0

(h/2),

Na ich podstawie wyznaczyć T

1

(h) i następnie obliczyć T

0

(h/4),

Na podstawie znanych wartości wyznaczyć T

1

(h/2), T

2

(h), itd.

17

W ten sposób na podstawie tablicy rzędu n otrzymuje się tablicę rzędu (n+1)

obliczając T

0

(h/2

n+1

) i na tej podstawie wyznaczając elementy T

k

(h/2

n+1-k

).

Obliczenia kończy się gdy różnica między T

n+1

(h) i T

n

(h) staje się mniejsza

niż założona wcześniej wartość.

MNM Całkowanie numeryczne M.Pyrz 04-2014

Metoda Romberga

( )

( )

0

0

0

0

0

2

3

4

1

1

1

1

2

3

2

2

2

2

2

2

2

h

h

h

h

T h

T

T

T

T

h

h

h

T h

T

T

T

 













 













 













 









 









 









⋯

⋯

18

Kryterium zatrzymania

( )

( )

1

n

n

T

h

T h

ε

+

−

≤

( )

( )

2

2

2

2

3

3

2

2

2

...

h

h

T

T

T h

h

T h

T

 





 





 





 

 

 

⋯

⋯

MNM Całkowanie numeryczne M.Pyrz 04-2014

Metody Gaussa

Całkowanie numeryczne

Idea: wybrać położenie N węzłów x

k

(i odpowiadających im współczynników

wagowych A

k

) w taki sposób, aby anulować błąd R

N

aproksymacji Î

N

.

We wzorze na całkowanie występuje 2N niewiadomych (x

k

i A

k

). Zatem

możemy spróbować je określić dla wielomianu interpolacyjnego stopnia

1

ˆ

( )

(

)

b

N

N

N

k

k

N

k

a

I

f x dx

I

R

A f x

R

=

=

≈

+

=

+

∑

∫

19

możemy spróbować je określić dla wielomianu interpolacyjnego stopnia

k=2N-1 (w którym występuje 2N współczynników).

Kwadratury metody Gaussa wykorzystywane są do obliczania całek w

granicach [-1,1], (-

∞

, +

∞

), [0,

∞

) …

W dalszej części rozważany będzie przypadek [a,b]=[-1,1] tj.

MNM Całkowanie numeryczne M.Pyrz 04-2014

1

1

( )

I

f x dx

+

−

=

∫

Metoda punktów Gaussa

Stosowana do całkowania wielomianów w granicach [-1,+1]
(duże zastosowanie w Metodzie Elementów Skończonych).
Całka obliczana jest jako suma wartości wielomianu wyznaczanych w

tzw. punktach całkowania, mnożonych przez współczynniki wagowe

1

1

( )

( )

n

i

i

i

I

f x dx

f

w

ξ

+

=

=

=

∑

∫

20

ξ

i

– położenie i-ego punktu Gaussa,

w

i

– współczynnik (waga),

n – liczba punktów całkowania.

Stosując

n

punktów całkowania można otrzymać dokładny wynik dla

wielomianu stopnia

2n-1

(lub niższego)

(na przykład do całkowania wielomianu stopnia 3,

potrzeba co najmniej 2 punktów Gaussa).

1

1

i

=

−

∑

∫

MNM Całkowanie numeryczne M.Pyrz 04-2014

Metoda punktów Gaussa

Położenie punktów całkowania (punktów Gaussa)

- symetryczne rozmieszczenie względem środka przedziału całkowania

- współczynniki (wagi) punktów położonych symetrycznie są jednakowe

Położenie

ξ

i

i wagi w

i

punktów Gaussa

21

Położenie

ξ

i

i wagi w

i

punktów Gaussa

Liczba punktów

Położenie punktów

Współczynniki wagowe

1

ξ

1

=0.0

w

1

=2.0

2

ξ

1

=

ξ

2

=

±

0.57735

w

1

=w

2

=1.0

3

ξ

1

=

ξ

3

=

±

0.774596

ξ

2

=0.0

w

1

=w

3

=0.555556

w

2

=0.888889

4

ξ

1

=

ξ

4

=

±

0.86113

ξ

2

=

ξ

3

=

±

0.33998

w

1

=w

4

=0.34785

w

2

=w

3

=0.65214

MNM Całkowanie numeryczne M.Pyrz 04-2014

Całkowanie wielomianu - przykład

Przypadek n=2. Wielomian stopnia 3 do scałkowania

Całkowanie dokładne prowadzi do wyniku

Jeżeli zastosujemy 2 punkty całkowania

ξ

1

=+k,

ξ

2

=-k , będziemy

stosowali wzór

f

a

a

a

a

( )

ξ

ξ

ξ

ξ

= +

+

+

1

2

3

2

4

3

I

f

a

a

=

=

+

−

+

z

( )

ξ

2

2

3

1

3

1

1

22

wynika stąd błąd obliczeń

Błąd zniknie (dla dowolnych wartości a

1

i a

3

) jeżeli spełnimy zależności:

czyli

2

1

3

(

)

( )

2 (

)

a

I

wf

k

wf k

w a

a k

=

− +

=

+

2

1

3

1

2 (1

)

2

(

)

3

a

I

I

a

w

a

wk

−

=

−

+

−

1

0

1

3

0

2

− =

−

=

w

wk

w

k

=

= ±

= ±

1 0

1

3

0 57735

.

.

MNM Całkowanie numeryczne M.Pyrz 04-2014

Całki wielokrotne

Analogiczne sformułowanie

Dobór punktów w zależności od stopnia

wielomianu względem poszczegόlnych

zmiennych

Carl Friedrich Gauss

(1777-1855)

1 1 1

p

n

m

+ + +

∑∑∑

∫ ∫ ∫

23

MNM Całkowanie numeryczne M.Pyrz 04-2014

(1777-1855)

1

1

1

1 1 1

( , , )

( ,

,

)

i

j

k

i

j

k

i

j

k

I

f x y z dxdydz

f

w w w

ξ η ζ

=

=

=

− − −

=

=

∑∑∑

∫ ∫ ∫