close

Вход

Забыли?

вход по аккаунту

код для вставкиСкачать
Richardson Extrapolation
ME6758, Dr. Ferri
Start with some integration rule, N(h), and some general trend of the error
as a function of h:
2
3
(1)
M  N (h )  K 1h  K 2 h  K 3 h  
Now, divide h by 2 (double the number of intervals)
h
h
h
M  N    K1   K 2  
2
2
2
2
3
h
 K3   
2
(2)
To eliminate the lowest order error term, take 2*Eq2 – Eq1


h
2 M  M   2 N    N ( h )   { K 1 h  K 1 h} 
2



 h2
 K 2 
2


 


 3
  K h 2    K  h   K h 3   
2
3
3

4 






 
 h2 
 h3 
 h

h



M   N    N    N (h)  K 2
 3K 3 



2
4 
2
 2





Call this N2(h)

h  h
h
N 2 ( h )  N     N    N ( h )   2 N    N ( h )
2  2
2

“weighted average”
Now,
 h2
M  N 2 (h)  K 2 
 2


 3
  3K  h   
3

4 



(3)
Halve h again (double the number of intervals):
3
h 1
2
3
M  N 2    K 2h 
K 3h  
32
2 8
(4)
Eliminate h2 terms by taking 4*Eq4 – Eq3
3
h
3
 3
4 M  M  3M  4 N 2    N 2 (h)   K 3  K 3  h  
4
2
8

M 
1
h 1
3
N 2    N 2 (h)  K 3 h  
3
8
2 3
4
Call this N3(h)
N 3 (h) 
Now,
h 1
h
N 2    N 2 (h)  N 2   
3
2 3
2
4
M  N 3 (h) 
1
8
K3 h
3
h
N 2    N 2 (h)
2
3
This is accurate
to order O(h3)

The O(h4) approximation is given by:
h
N 4 (h)  N 3   
2
h
N 3    N 3 (h)
2

7
h 1
N 3    N 3 (h)
7
2 7
8
The O(h5) approximation is given by:
h
N 5 (h)  N 4   
2
In general:
h
N 4    N 4 (h)
2
h
N j ( h )  N j 1   
2
15

1
h
N4  
N 4 (h)
15
 2  15
16
h
N j 1    N j 1 ( h )
2
2
j 1
1
This is accurate
to order O(hj)
Construct a Table
O(h)
O(h2)
O(h3)
O(h4)
N1(h)
N1(h/2)
N2(h)
N1(h/4)
N2(h/2)
N3(h)
N1(h/8)
N2(h/4)
N3(h/2)
N4(h)
Romberg Integration
Start with composite trapezoidal rule with m-intervals, (m+1) points
m 1


ba 2
M   f ( x ) dx   f ( a )  f ( b )  2  f ( x j )  
h f  (  )
a
2
12

j 1


b
h
h = (b-a)/m
xj =a + j*h
a<<b
Perform n trapezoidal integrations using m1 = 1 interval, m2 = 2 intervals,
m3 = 4 intervals, … mn = 2n-1 intervals. In each case hk = (b-a)/mk = (b-a)/ 2k-1
k 1


2
1
ba 2
 f ( a )  f ( b )  2  f ( a  ih k )  
M  
f ( x ) dx 
h k f  (  k )
a
2 
12

i 1


b
hk
Call this Rk,1
R1,1
R 1,1 
h3
h2
h1
R2,1
h1
[ f ( a )  f ( b )] 
2
R 2 ,1 
h2
R 2 ,1 
h2
2
2
ba
R3,1
[ f ( a )  f ( b )]
2
[ f ( a )  f ( b )  2 * f ( a  h 2 )]
[ f ( a )  f ( b )]  [ h 2 f ( a  h 2 )] 
h1/2
just R1,1/2
1
2
[ R 1,1  h1 f ( a  h 2 )]
h3
h2
h1
R1,1
R2,1
R 3 ,1 
h3
2
R3,1
[ f ( a )  f ( b )  2 f ( a  2 h 3 )]  h 3 [ f ( a  h 3 )  f ( a  3 h 3 )]
just R2,1/2
R 3 ,1 
1
2
R 2 ,1  h 2 [ f ( a  h3 ) 
h2/2
f ( a  3 h 3 )]
k 2

2
1
R k ,1   R k  1,1  h k  1  f ( a  ( 2 i  1) h k
2
i 1




)

odd multiples
of hk
R1,1
new point
R2,1
new point
new point
R3,1
h4
R4,1
h5
R5,1
R6,1
Rk,1
2 endpoints, 2k-1-1 interior points, 2k-2 additional points
odd multiples of hk

0
Example
k  1, 2
k 1
k  2, 2
k  3, 2
k  4, 2
k  5, 2
 1, h1   , R1,1 
k 1
k 1
 2 , h 2   / 2 , R 2 ,1 
 8 , h 4   / 8 , R 4 ,1 
1
 16 , h 5   / 16 , R 5 ,1 
1
 sin(
5
16
k  6, 2
k 1
[sin( 0 )  sin(  )]  0
2
 4 , h 3   / 4 , R 3 ,1 
k 1
k 1

)  sin(
7
16
sin x dx  2
)  sin(
9
16
2
2
1
2
1
2
[ R1,1   sin(
[ R 2 ,1 
[ R 3 ,1 
[ R 4 ,1 
)  sin(


[sin(
2
[sin(
4

8
11 
16


2
)  sin(
3
4


)]]  1 . 8962
4
)  sin(
3
)  sin(
8
)  sin(
16
)  sin(

2
8
[sin(
)] 
3
5
)  sin(
8
)
16
13 
)  sin(
16
 32 , h 6   / 32 , R 6 ,1  1 . 9984
Converges, but rather slowly…
15 
16
)]]  1 . 9936
7
8
)]]  1 . 9742

b
 a f ( x ) dx  R k ,1 
 K i hk
2i
i 1


b
 a f ( x ) dx  R k  1,1 

i 1
2
 K 1hk 

2i
K i hk 1 
 K i hk
2i
i2
2i
 Ki
i 1
hk
2
2i

2

K 1hk

4
2i
 Ki
i2
hk
4
i
Subtract top equation from 4 times bottom equation. After algebra, get
b
a
f ( x ) dx 
4 R k  1,1  R k ,1
3

K i  1  4 i 1


3  4 i 1
i2


 h 2i
 k

Re-write as
R k  1,1  R k ,1   K i

 
 a f ( x ) dx   R k  1,1 
3

 i2 3
b
 1  4 i 1

 4 i 1


 h 2i
 k

Now, extrapolate again to eliminate O(hk4) terms to obtain an O(hk6) result, etc
Define:
Generalize to:
R k , 2  R k ,1 
R k ,1  R k  1,1
3
R k , j  R k , j 1 
R k , j  1  R k  1, j  1
4
Obtain from
a single
composite
trapezoidal
integration
R1,1
j 1
Error:
O(hk2j)
1
Construct Romberg Table
R2,1
R2,2
R3,1
R3,2
R3,3
R4,1
R4,2
R4,3
R4,4
Rn,1
Rn,2
Rn,3
Rn,4
I=Rn,n + O(hn2n)
…
Rn,n
Obtain by simple “averaging”

0
Example
sin x dx  2
0.00000000
1.57079633
2.09439511
1.89611890
2.00455976
1.99857073
1.97423160
2.00026917
1.99998313
2.00000555
1.99357034
2.00001659
1.99999975
2.00000001
1.99999999
1.99839336
2.00000103
2.00000000
2.00000000
2.00000000
Error in R6,6 is only 6.61026789e-011 % !!!!
h6 = /32= 9.817e-002,
h612 = 8.017e-013
2.00000000
Accurate to
O(h612)
1/--страниц
Пожаловаться на содержимое документа