• Back to MAIN


    • Let's define J(f;a;b) as the INTEGRAL of the function f over the closed interval T=[a;b], we look here for
    a J* holding abs(J-J*) < e for any given real e > 0

    So..if f = x4 + x3 .. then J(f;0;2) = 10.4. In case we have a finite number of f-values {f(xi)} where xi in [a;b],
    then J* must be of the form:

    J* = SUM wif(xi), where the wi are the weights of the quadrature.
    Decomposing the interval regularly into n bits.. we obtain easily the relation xi = a + ib ...where i in (0,...,n), and let us
    further define xr.5 as the center of the interval [xr;xr+1]. The interval definitions allow us to define the trapezoid-rule TA

    TA(h) = h*[0.5*f(x0)+SUM[f(xi);0;n-1]+0.5*f(xn)]

    and the center-rule MA

    MA(h) = h*SUM[f(xi.5);0;n-1]

    As T(h/2) = (by doubling n) = 0.5*h*[0.5*f(x0)+f(x0.5)+..+f(x(n-1).5)+0.5*f(xn)]=
    0.5*h*[0.5*f(x0)+f(x1)+..+0.5*f(xn)]+0.5*h*[f(x0.5)+..+f(x(n-1).5]=0.5*[MA(h)+TA(h)]..which directly implies that

    T(h/2) = 0.5[TA(h)+MA(h)]

    ..so, we see here that doubling of the precision implies doubling of the functional evaluations.., intuition tells us here already,
    that f is not allowed to be "too wild"..eg should not have kind of "asymptotical misbehaviour"...,but even then, the
    trapezoid and center - rules converge quite slowly.

    A trivial solution in k4 representing the trapezoid rule might be :

    mod:{x-y xbar x}
    xexp:{exp y*log x}
    xe3:{xexp[x;3f]}
    rot:{,/|(0;mod[x;#y])_y}

    Ta:{[f;a;b;n]h*+/(2 _ r),0.5*2 # r:rot[-1;f (!n+1)*h:(b-a)%n]}

    So..using our f = x4 + x3 we get : Ta . ({(x+1)*xe3 x};0;2;20) returns : 10.43666 ..and Ta . ({(x+1)*xe3 x};0;2;20) returns: 10.4
    But this was all just scratching a bit on the surface of this problem..

    Leonhard Euler (1707-1783) and Colin MacLaurin (1698-1746) independently investigated the error of the trapezoid-rule in the
    asymptotic case. Let's look at the general MacLaurin serie..

    f(a + h) = f(a) + h f'(a)/1! + h2f''(n)/2!+ ... + hnf(n)/n! + ..
    So, in case this representation of f exists in [a;b], we deduct easily:

    f(a + h) - f(a) = [1 + h f'(a)/!1 + h2f''(n)/!2 + ... + hnf(n)/!n] - 1 = (ehD - 1)f(x)
    where D is the differential d/dx, and consequently

    f(x) = (ehD - 1)-1[f(a + h) - f(a)]=SUM[Bk(hD)k-1/!k;k=0;k=infinity] ,
    where Bk are the Bernoulli-numbers. Jacob Bernoulli (1654-1705)

    A friendly hint:..we can obtain the Bernoulli-numbers by solving the infinite equation-system:

    x0 = 1
    2x1 + 1 = 0
    3x2+3x1 + 1 = 0
    4x3+6x2+4x1 + 1 = 0
    etc..

    ..which k4 blows away by the one-liner

    {|(x#-1f)$+!:1f*(-x)#/:(x#0),/:1 _ {1 _ -1 _ x}'x {({+':x}@x),1}\ 1 1}

    this "lazy" way works until n = 28 (numerical extinction process of the matrix-inversion spoils it, don't forget we handle here
    an unfriendly non-unitary matrix). BUT:.. a K recursion is stable for ALL n ..(this K-recursion should be now
    a fair exercise for the non-K beginner reader...it is even shorter than our "lazy" way).. Lets "accept" for a while our
    simple solution..returning us(..and beginning with x1):

    (-0.5,0.1666667,0,-0.03333333,0,0.02380952,0,-0.03333333,0,0.07575758,0,-0.253113,0 ... etc.)

    Analytically all B2n+1 are 0 for all n>0.

    From the equation for f(x) above follows(if we start the sum with index k=2), that:
    f(x) = {(hD)-1 - 0.5 + SUM[Bk(hD)k-1/!k;k=2;k=infinity]}(f(x + h) - f(x)) .. solving this for f(x + h) + f(x)
    and because D-1 is an integral, we will get:

    0.5*h*[f(x) + f(x + h)] = T(h) = J[f;x;x + h] + SUM[ [B2kh2k;k=1;k=infinity](f2k-1(x + h) - f2k-1(x)/(2k)!);k=1;k=infinity) ],
    as this works for all h in [a;b] we finally will obtain the socalled Euler-MacLaurin sum-formula:
    T(h) = J[f;a;b] + SUM[ukh2k;1;N] + RN+1(h) .. where uk = B2kh2k(f2k-1(b) - f2k-1(a))/(2k)!) and RN+1 a residual term.
    Here we see already, that in the general case where uk is not equal to 0 (asymptotic behavior), the error of a pure
    trapezoid method is in the order of h2, and because we have here only even exponents, we can extrapolate the trapezoid
    sums by the Romberg-Schema as a special case of the algorithm of Watson Neville (1886-1965):

    Tm+k,k = (4kTm+k,k-1 - Tm+k-1,k-1)/(4k - 1) for all m > 0, k > 1
    and for all Tm,0 , where m > 0, we identify as our old trapezoid rule


    leading us to the triangular schema:

    T+1,0
     T2,0 , T+2,1
     T3,0 , T3,1 , T+3,2
     T4,0 , T4,1 , T4,2 , T+4,3
     etc..

    where the T+ sequence represents the convergence-path. One possible answer in k4 to this problem can be
    (using a locale romberg.k):

    \d .romberg
    /assuming that a,b,f and n are given

    xe4:.r.xe4
    p2:.r.p2
    p4:.r.p4

    tr:{x+(x-y)%-1+p4 z}
    ts:{[f;a;b;n]l:0.5*last u:f@a+(h:(b-a)%r)*i:!1+r:p2 n;h*+/@[u;(*i),last i;:;l,0.5**u]}

    /romberg schema
    s:{ts[f;a;b;x]}'h2:!v:1+n;
    m:v#'v#0f;m[;0]:s;i:0;do[n;k:i+1;j:0;do[k;v:tr . (m[k,i;j],j+1);m[k;j+1]:v;j+:1;];i:i+1;]

    /end of locale

    \c 23 155
    \d .r

    expo:{exp y*log x}
    xe4:{expo[4f;x]}
    p2:31 {x*2}\1
    p4:15 {x*4}\1

    q)setx:{.["S"$".romberg.",string x;();:;y]}
    \P 11
    romberg:{[f;a;b;n]`f`a`b`n setx' (f;a;b;n);.:"\\l romberg.k";.romberg.m}

    ..where i is the depth, and where the function romberg is returning the convergence path. The Romberg schema terminates
    if we find an [i] where the absolute difference of 2 neighbor-diagonal elements are less than a given epsilon, which is a
    trivial K exercise.

    So, in case of g1 = x-1 integrated over the interval [1;2], last .r.romberg . ({%x};1f;2f;4)
    returns (0.69339120221 0.69314765282 0.6931471943 0.69314718307 0.69314718192

    In case of g2 = 4(x2 + 1)-1 over [0;1] .. last .r.romberg . ({4f%1+x*x};0f;1f;5)
    returns (3.1414298932 3.1415926536 3.1415926537 3.1415926536 3.1415926536 3.1415926536)
    In case of g3 = x10e4x3-3x4 over [0;2] .. last .r.romberg . ({expo[x;10f]*exp x*x*x*4-3*x};0f;2f;7)
    returns (7.2583950704 7.2583951681 7.2583951704 7.2583951705 7.2583951633 7.2583951616 7.2583951612 7.2583951611)

    In case of g4 = 1 + e-x sin 8x2/3 over [0;2] .. last .r.romberg . ({1+(exp (-x))*sin 8*expo[x;2%3]};0f;2f;8)
    returns (2.0159008916 2.0161744346 2.0161897633 2.0161928752 2.0161936173 2.0161938008 2.0161938465 2.0161938579 2.0161938608)

    In case of g5 = (14x-11x2)e-2x over [0;2] .. last .r.romberg . ({(x*14f-11f*x)*exp -2f*x};0f;2f;5))
    returns (1.0797214333 1.0842349364 1.0842600515 1.084260394 1.0842604075 1.0842604089)

    In case of gauss = (2*Pi)-0.5e-0.5x2 over [-3;3] .. last .r.romberg . ({(%sqrt 2*acos -1f)*exp -0.5*x*x};-3f;3f;7)
    returns (0.99729533603 0.99730019966 0.99730020393 0.99730020394 0.99730020394 0.99730020394 0.99730020394 0.99730020394

    The last values are the correct values, we see, the Romberg schema converges quickly, usually we for i=5 or i=6 we
    obtain already the expected result, further iterations can distort the result by round-off errors.

    Basically, the Romberg schema eliminates the first occurrence of the error u1h2 in the Euler Mac-Laurin sum formula.
    The Romberg algorithm eliminates successively all error terms h4,h6..so that Tm+k,k converges to J[f;a;b] plus an error
    term of the order h2k+2.

    The Romberg schema relies on the fact, that the derivatives are finite in the end-points (without explicitly calculating them).
    The Romberg schema will fail in case of an integrand having a divergent slope at one or both boundary ends.

    One such an example would be: h(x) = (1 - x2)0.5 over the interval [0;1]

    But we should not waste to much time with investigations of Romberg's algorithm regarding those singularities. There are far more
    suited methods, for that particular case. We will not discuss these methods in this chapter.

    As we saw, the Romberg schema is a fairly simple and reliable method, but can turn out to be quite CPU-expensive, if we really need a
    very precise result.
    Setting hj-1 = h and hj = 0.5h, we obtain from Tj,1 = (4Tj,0 - Tj-1,0)/3 and from the equation above that T(h/2) = 0.5[TA(h)+MA(h)]
    follows Simpson's formula,Thomas Simpson (1710-1761):

    Tj,1 = [4T(h/2) - T(h)]/3 = [T(h) + 2M(h)]/3
    remembering that the Tj,1 approximate the integral with an error of the order h4

    for a general h = (b - a)/n S(h) = h[ 0.5(f(x0) + f(xn)) + SUM[f(xj);j=1;n] + 2 SUM[f(xj.5);j=0;j=n-1] ]/3


    The applied methods so far did not take into consideration some aspects like the curvature, this means the romberg
    algorithm can be a very uneconomic waste regarding evaluations.

    This should improve with an adaptive method, which means that we take (one possible approach) as the 1st approximation the
    trapezoid value, and for the 2nd approximation we take the Simpson value: so, for hj:=bj - aj and mj=0.5(aj+bj) we set
    Tj=0.5hj[f(aj) + f(bj)] and Sj=[Tj + 2hjf(mj)]/3, substituting Tj, we finally get:

    Sj = h[f(aj) + 4jf(mj) + f(bj)]/6
    This formula is reminding us of something else, it is a method originally stated from the famous astronomer
    Johannes Kepler (1571-1630). Truth is, Simpson basically re-used Kepler's method to calculate the integral.
    Applying Kepler's rule on the interval [xi,xi+2] we will get for even n, starting with 0 and ending with n-2, and h = 2h*

    Sj = h*[f(xj) + 4f(xj+1) + f(xj+2)]/3
    j in [0;n-2]


    ..which is the more common form.

    Resolving that in k4 can look like this:

    expo:{exp y*log x}
    xe2:{expo[2f;x]}

    simp:{[f;a;b;n](e%3*r)*+/+/(1,4,1)*/:(f a+(! _ r+1)*(e:b-a)%r)@0 1 2+/:2*!_ .5*r:xe2 n}
    Let's try a nastier integral of f = sin (1+x3)0.5 + e-3x + e-40(x - 0.5)2 over [0;1] .. simp . (f;0;1;13)
    where f:{(sin sqrt 1+x*x*x)+(exp x*-3f)+exp -40f*r*r:x-0.5} returns : 1.486643620753188

    In fact, this first method does not really beat our Romberg schema .. which is not much surprising, as the Simpson
    method basically replaces the linear approach by a parabola approach. Our sample function was obviously not very
    appropriate for this case, too:(the x-axis is in 0.0001 units). On the other hand, the simp above does too many function
    evaluations.(as one weakness of that proposed program solution)

    7

    so lets look at other and better methods..

    One other method are the Newton-Cotes, is using the approach of a piecewise weighting of n-degree polynomials.
    This method relies further on a unique representation of the integral J[f;xi;xi+j-1] by the sum :

    t-1(xi+j-1 - xi)SUM[rsyi+s;s = 0;s = n - 1] + errj(f)

    and where xs = xi + h*s,(s = 0,..., j - 1), and h = (j - 1)-1(xi+j-1 - xi), and where yi = f(xi).
    In standard math books we find all those values listed..here and excerpt going up to the polynomial degree j = 11.
    We talk here about the closed Newton - Cotes.

    j h-factor r0 r1 r2 r3 r4 r5 r6 r7 r8 r9 r10
    2 1 1 1
    3 1/3 1 4 1
    4 3/8 1 3 3 1
    5 2/45 7 32 12 32 7
    6 5/288 19 75 50 50 75 19
    7 1/140 41 216 27 272 27 216 41
    8 7/17280 751 3577 1323 2989 2989 1323 3577 751
    9 4/14175 989 5888 -928 10496 -4540 10496 -928 5888 989
    10 9/89600 2587 15741 1080 19344 5788 5788 19344 1080 15741 2587
    11 5/299376 16067 106300 -48525 272400 -260550 427368 -260550 272400 -48525 106300 16067


    Friendly hint: All Newton-Cotes for even j = 2q have a worse precision than the method using the degree of 2q-1.

    For j = 2 we find our Simpson rule again. All these coefficients can be calculated by using Lagrange-polynomials.
    For training reasons we concentrate here on the fix polynomial approximation of degree 7.

    We take here as sample an integral from physics, expressing the heat capacity under a constant volume

    V:Cv = 2.25R rm3 = J[x4 sinh -2 0.5x;0;rm]

    where R is the universal gas constant, T = temperature, rm = dbT-1 and where db is the debye temperature of a solid
    (the Debye temperature is the temperature of a crystal's highest normal mode of vibration, it is db = hvmk-1
    h is Planck's constant, k is Boltzmann's constant, and vm is the Debye frequency.)
    So f expressed in k4 : {r*r:(x*x)%sinh x%2f}, where sinh:{0.5*(%r)-r:exp@-x}. f converges to 0
    if x goes against 0, x-unit = 0.01 in picture below

    7

    We search for the integral value of J[f;0;30] .. using the 7th grade polynomial.
    The calculation using simp brings us by : simp . (f;1.0e-14;30;9) to the value: 103.903030088336

    One possibility in k4 for such a 7th grade approximation would be:

    w:41 216 27 272 27 216 41;hq:%140
    g7:{[f;a;b;q;n]!m:_ .5+((qn:expo[q;n*1f])-1f)%q - 1;cf:(-!m)+'(m,q)#!q*m;+/+/hq*h*w*/: ..
    .. (f s:a,((1 _ ! _ qn)*h:(b-a)%qn),b*1f)@cf}

    g7 . (f;1.0e-14;30;7;5) would give us 103.903030087795

    the graphical integration behavior of g7:..(in g7 we just replace the return: +\+/+ hq*h*w*/:(f s .. etc

    7

    g7 is a pure demo resp. school version, we see many reasons why we would not use g7 in a production system.

  • Back to MAIN


    • ©++ MILAN ONDRUS