Derivations for the moving
parabolic approximation (in Maple)




This is a copy of an interactive session with Maple where we compute the main formulas used for the parabolic approximations.

    |\^/|     Maple V Release 4 (ETH Zuerich)
._|\|   |/|_. Copyright (c) 1981-1996 by Waterloo Maple Inc. All rights
 \  MAPLE  /  reserved. Maple and Maple V are registered trademarks of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> S2 := sum( ( x[i-j] - (a0[i]+a1[i]*j+a2[i]*j^2) ) ^ 2, j=0..m-1 );
                   m - 1
                   -----
                    \                                         2 2
             S2 :=   )   (x[i - j] - a0[i] - a1[i] j - a2[i] j )
                    /
                   -----
                   j = 0

> vars := {a0[i], a1[i], a2[i]}:
> sol := solve( map2(diff,S2,vars), vars ):

The approximation problem is posed as a least squares problem. Here the sum of squares is stored in S2 and its derivatives with respect to the unknown variables (vars) are equated to 0. The solution contains three different explicit summations. It is useful to represent them with symbols. At this point the solution is too large and we do not show it. It needs to be simplified. We build a pattern to do this substitution and to keep as a comment.

> subspat := [sum(x[i-j],j=0..m-1) = sx,
>             sum(j*x[i-j],j=0..m-1) = sjx,
>             sum(j^2*x[i-j],j=0..m-1) = sj2x ];
subspat :=
     m - 1                m - 1                   m - 1
     -----                -----                   -----
      \                    \                       \     2
    [  )   x[i - j] = sx,   )   j x[i - j] = sjx,   )   j  x[i - j] = sj2x]
      /                    /                       /
     -----                -----                   -----
     j = 0                j = 0                   j = 0

> sol := factor( subs( subspat, sol ));
                        2
                  3 sx m  - 12 sjx m - 3 sx m + 10 sj2x + 2 sx + 6 sjx
sol := {a0[i] = 3 ----------------------------------------------------, a1[i]
                                   m (m + 2) (m + 1)

                 3           2          2
     = -6 (6 sx m  - 32 sjx m  - 21 sx m  + 60 sjx m + 30 sj2x m + 21 sx m

     - 22 sjx - 30 sj2x - 6 sx)/(m (m - 1) (m - 2) (m + 2) (m + 1)),

                   2
               sx m  - 3 sx m - 6 sjx m + 6 sj2x + 2 sx + 6 sjx
    a2[i] = 30 ------------------------------------------------}
                      m (m - 1) (m - 2) (m + 2) (m + 1)

Factoring the solutions provides some significant simplification. There are many common subexpressions, so it is obviously useful to optimize these computations.

> readlib(optimize):
> Comp_as := [optimize( [op(sol)], tryhard )];
                               2                             1
Comp_as := [t15 = sjx m, t8 = m , t14 = sx t8, t13 = -----------------,
                                                     (m + 1) (m + 2) m

                                             t13
    t12 = 6 sjx + (-3 m + 2) sx, t11 = ---------------,
                                       (m - 2) (m - 1)

    a0[i] = 3 (3 t14 - 12 t15 + 10 sj2x + t12) t13,

    a1[i] = -6 ((30 m - 30) sj2x + (-32 t8 - 22 + 60 m) sjx

     + (-21 t8 - 6 + (6 t8 + 21) m) sx) t11,

    a2[i] = 30 (t14 - 6 t15 + 6 sj2x + t12) t11]

The optimized code now contains assignments to several temporary variables (t15, t8, ... etc.) To produce the OpenMath block to do these computations we need to collect all new temporary variables and produce a Declaration statement. Also, equations should be converted to Assign statements.

> Locs := {op(map2(op,1,Comp_as))} minus vars;
                     Locs := {t8, t11, t12, t13, t14, t15}

> Compute_as := Block( Declaration( [op(Locs)] = [Local,Float] ),
>       seq( Assign( op(1,z), op(2,z)), z=Comp_as ));
Compute_as := Block(

    Declaration([t8, t11, t12, t13, t14, t15] = [Local, Float]),

                                    2
    Assign(t15, sjx m), Assign(t8, m ), Assign(t14, sx t8),

                        1
    Assign(t13, -----------------),
                (m + 1) (m + 2) m

    Assign(t12, 6 sjx + (-3 m + 2) sx),

                      t13
    Assign(t11, ---------------),
                (m - 2) (m - 1)

    Assign(a0[i], 3 (3 t14 - 12 t15 + 10 sj2x + t12) t13),

    Assign(a1[i], -6 ((30 m - 30) sj2x + (-32 t8 - 22 + 60 m) sjx

     + (-21 t8 - 6 + (6 t8 + 21) m) sx) t11),

    Assign(a2[i], 30 (t14 - 6 t15 + 6 sj2x + t12) t11))

The block is now ready to be inserted in the main function.