Re: Romberg Interpolation
Soroosh Yazdani
be just the values you get using Trapazoid approximation.
Richardson extrapolation of the Trapazoid rule should always converge to
the actual integral, for integrable functions. Romberg integration
is supposed to be applying Richardson extrapolation to Trapazoid rule,
however, the books implementation might have a bit of twist that
prevents the method from worknig all the time. I don't have the book on me
right now, but I will look into it.
Soroosh
On Mon, Oct 30, 2006 at 10:36:51PM -0800, mvish...@berkeley.edu wrote:
> Hey Soroosh,
>
> Thanks for you help in office hours about Romberg interpolation. I was
> wondering if there some property of the function f that must be satisfied
> for the Romberg interpolation on f to converge to the integral value? For
> some reason, as I reduce my n, my solution gets closer to the answer
> (which would defy the purpose of iteratively "refining" the solution). I
> wrote my own code (instead of the algorithm in the book), and I'm pretty
> sure that it is correct.
>
> If you don't have to look at this, that is fine, because I know you're
> busy too. I spent over two hours thinking about why this wasn't working,
> and so I abandoned it and have three methods (Composite Simpsons,
> Trapezoidal, and Midpoint). Out of curiosity, I was wondering why this
> isn't behaving.
>
> Thanks a lot!!!
>
> Meghana
>
>
> Here's the code:
>
> function [R, efficiency, diff] = Romberg(f, a, b, n);
>
> f=fcnchk(f);
> R=zeros(n);
> h=(b-a);
> exact=4/3;
> sumo=0;
>
> R(1,1)=h*(f(a)+f(b))/2;
>
> for k=2:n,
> for i=1:2^(k-2)
> m = a + ((2*i-1)*(h/2^(k-1)));
> sumo = sumo + f(m);
> end
> R(k,1) = (R(k-1,1) + (h/2^(k-2))*sumo)/2;
> end
>
> for k=2:n,
> for j=2:k,
> R(k,j) = R(k, j-1) + ((R(k,j-1)-R(k-1,j-1))/(4^(j-1)-1));
> end
> end
>
> diff = abs(R(n,n)-exact);
> efficiency=n;
>
>
> ----------------
>
> >> [a,b,c]=Romberg('sqrt(abs(x))', -1, 1, 5)
>
> a =
>
> 2.0000 0 0 0 0
> 1.0000 0.6667 0 0 0
> 1.2071 1.2761 1.3168 0 0
> 1.6401 1.7845 1.8183 1.8263 0
> 2.0113 2.1351 2.1584 2.1638 2.1651
>
>
> b =
>
> 5
>
>
> c =
>
> 0.8318
>
>
> ----------------
>
> >> [a,b,c]=Romberg('sqrt(abs(x))', -1, 1, 10)
>
> a(n, n) = 2.6387 (I cut out the matrix to save space).
>
>
> b =
>
> 10
>
>
> c =
>
> 1.3054
>
>
>
>
>