Adaptive rational interpolation for point values

Journal of Computational and Applied Mathematics (

)

–

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Adaptive rational interpolation for point values✩ Francesc Aràndiga Departament de Matemàtica Aplicada, Universitat de València, Spain

article

info

Article history: Received 16 March 2018 Received in revised form 18 August 2018 Keywords: Multiscale decomposition Interpolation Adaptive rational interpolation Non linear weights

a b s t r a c t G. Ramponi et al. introduced in Carrato et al. (1997,1998), Castagno and Ramponi (1996) and Ramponi (1995) a non linear rational interpolator of order two. In this paper we extend this result to get order four. We observe the Gibbs phenomenon that is obtained near discontinuities with its weights. With the weights we propose we obtain approximations of order four in smooth regions and three near discontinuities. We also introduce a rational nonlinear extrapolation which is also of order four in the smooth region of the given function. In the experiments we calculate numerically approximation orders for the different methods described in this paper and see that they coincide with those that have been obtained theoretically. We also present reconstructions in 1d and 2d with which we reach the same conclusions. © 2018 Elsevier B.V. All rights reserved.

1. Introduction Image zooming can be performed using the tensor product of one dimensional interpolatory techniques. So we can consider an image of n × n pixels as a function defined on [0, 1] × [0, 1], and the value of the pixel (i, j) represents f (xi , yj ), where {xi , yj }1≤i,j≤n is a discretization of the square [0, 1] × [0, 1]. Linear techniques, data independent, as the Lagrange polynomial interpolation, lead to fast and efficient algorithms. However, the quality of the zoomed image may not be satisfactory in the regions where the stencils used to construct the interpolatory function contain a singularity of the function. Increasing the order of the interpolatory technique does not solve the problem. In this case the number of stencils that cross the singularity grows, increasing the region of poor accuracy. The use of non-linear interpolatory techniques, data-dependent, produces an improvement in the visual quality of the zoomed image. The interpolatory technique ENO (Essentially Non-Oscillatory) introduced by A. Harten et al. [1] performs a previous choice of the stencil, so that the interpolant is constructed, if it is possible, taking information from regions where the function is smooth. The WENO (Weighted ENO) interpolatory technique introduced by Liu et al. [2] is an improvement of the ENO technique. This technique builds the interpolant using a convex combination of all the approximations computed with stencils containing the interval where we want to obtain an approximation. The key lies in a wise assignment of weights to the convex combination. These weights must be chosen so that polynomials corresponding to singularity-crossing stencils should have a negligible contribution to the convex combination. The rational interpolation introduced by S. Carrato and G. Ramponi [3–6] is a second order non-linear interpolatory technique where the convex combination is given by the end-points of the interval where we want to obtain an approximation. In this paper we describe an improvement on this nonlinear/adaptive rational interpolator. We analyze this method and we see that it has order of approximation two. We define a related procedure with new weights inspired in the WENO ✩ This research was partially supported by Spanish MCINN, Spain MTM2014-54388-P and MTM2017-83942-P. E-mail address: [email protected]. https://doi.org/10.1016/j.cam.2018.08.052 0377-0427/© 2018 Elsevier B.V. All rights reserved.

Please cite this article in press as: F. Aràndiga, Adaptive rational interpolation for point values, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.052.

2

F. Aràndiga / Journal of Computational and Applied Mathematics (

)

–

weights and obtain a reconstruction that is four order accurate. Our proposal leads to a reconstruction that fulfills the WENO properties and therefore our method it is comparable to the WENO interpolation. In fact, adaptivity comes from the fact that the weights are chosen to meet the WENO properties. The importance of using an efficient interpolatory technique is closely related to image compression. Harten’s multiresolution framework [7,8] uses a prediction operator which goes from a lower to a higher resolution level. The prediction operator involves the application of an interpolatory scheme to the lower resolution level. So, the use of an efficient interpolatory technique as prediction operator leads to larger compression capabilities. This paper is organized as follows. In Section 2 we describe the rational interpolation introduced by S. Carrato and G. Ramponi [3–6]. We also analyze its order of approximation. In Section 3 we extend this technique to obtain order four, analyzing the problems that occur close to jumps of the function. We propose new weights in Section 4. We prove that with them, we obtain reconstructions of order three close to the discontinuities of the function. In Section 5 we present rational extrapolators which are also of order four in smooth areas and of order 2 close to jumps of the function. In Section 6 we calculate numerically approximation orders for the different methods presented and see that they coincide with those that have been obtained theoretically in the previous sections. We also present reconstructions in 1d and 2d. 2. Rational interpolation of order two N

Let be X k = {xki }i=k0 a uniform discretization of the interval [0, 1], where k specifies the level of resolution, i.e. xki = i · hk , N hk = 1/Nk . Let us assume that the data we have, {fik }i=k0 , represents the point values of a function f (x) in X k , i.e. fik = f (xki ). The problem of point-value interpolation consists in approximating the point values of f (x) in a higher resolution level. N 1 k+1 = i · hk+1 , i = 0, . . . , Nk+1 , Nk+1 = 2Nk . Note that in this case We consider a finer grid, X k+1 = {xki +1 }i=k+ 0 , xi 1 k+1 k k x2i−1 = (xi−1 + xi )/2. The rational interpolation proposed in [3] consists in approximating f2ik+ −1 by a weighted average k k between fi−1 and fi , i.e. 1 k k fˆ2ik+ −1 = ωi−1 fi−1 + ωi fi , with ωi−1 + ωi = 1.

(1)

instead of the linear approximation 1 fˆ2ik+ −1 =

1 2

k fi− 1 +

1 2

1 2 fik (= f2ik+ −1 + O(hk )).

(2)

The weights proposed in [3] are the following

ωi0−1 = ωi0

=

k k 2 1 + α (fi− 1 − fi+1 ) k k 2 k k 2 2 + α (fi− 2 − fi ) + (fi−1 − fi+1 )

(

k k 2 1 + α (fi− 2 − fi ) k k 2 2 + α (fi−2 − fik )2 + (fi− 1 − fi+1 )

(

k

), ),

(3)

and k k k k 2 2 1 + α (fi− 1 − fi+1 ) + (fi − fi+1 )

(

ω

1 i−1

ωi1

= =

2+α 2+α

(∑1

k s=0 ((fi−1 k (fi− 1 1 k ((f i−1 s=0

1+α (∑

(

)

), k k k 2 2 − fi+ 1−3s ) + (fi − fi+1−3s ) ) ) k k k 2 2 − fi− 2 ) + (fi − fi−2 ) ), k k k − fi+1−3s )2 + (fi − fi+1−3s )2 )

(4)

k where α is a variable parameter. Note that if we take α = 0, the approximation is, in the two cases, the average between fi− 1 k and fi . In [3] the authors propose to use α = 0.01.

2.1. Order of accuracy 1 r Let us examine the order of accuracy that we get when we approximate f2ik+ −1 by the expression given in (1). Consider qi (x) k k the polynomial of degree r interpolating f (x) in the stencil {xi , . . . , xi+r }, i.e. k qri (xki+l ) = fi+ l,

l = 0, . . . , r ,

and let be pri = qri (xk2i+−11 ). Since p1i−1 = obtain

(5) 1 0 p 2 i−1

1 l 2 + 12 p0i = f2ik+ −1 + O(hk ), if we prove ωi−s =

1 k+1 l 0 l 0 2 fˆ2ik+ −1 = ωi−1 pi−1 + ωi pi = f2i−1 + O(hk ).

1 2

+ O(h2k ), l, s = 0,1, we hope to (6)

To see this, first we need the following proposition. Please cite this article in press as: F. Aràndiga, Adaptive rational interpolation for point values, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.052.

F. Aràndiga / Journal of Computational and Applied Mathematics (

)

–

3

Proposition 1. Let us assume that f is smooth in [xki−2 , xki+1 ]. Then

ωil−s =

1 2

+ O(h2k ), l, s = 0, 1.

(7)

Proof. If f is smooth, a Taylor expansion leads to: f (xki+n ) − f (xki+m ) = (n − m)f ′ (xki )hk + O(h2k ) = O(hk ), n ̸ = m.

(8)

In this way, if we use (8) in (3) and (4), we obtain:

ωil−s =

1 + α O(h2k ) 2 + α O(h2k )

, l, s = 0, 1

(9)

and

ωil−s −

1 2

=

1 + α O(h2k ) 2+α

−

O(h2k )

1 2

2α O(h2k ) − α O(h2k )

=

2(2 + α O(h2k ))

= O(h2k ), l, s = 0, 1.

□

Proposition 2. If f is smooth in [xki−2 , xki+1 ] then k+1 1 l 0 l 0 2 fˆ2ik+ −1 = ωi−1 pi−1 + ωi pi = f2i−1 + O(hk ),

l = 0, 1.

(10)

Proof. We have 1 ∑

1 ∑

ωil−s p0i−s − p1i−1 =

s=0

=

1 ∑

=

ωil−s p0i−s −

1 ∑ 1

2

s=0

ωil−s p0i−s − f2i−1 − ωil−s −

s=0

1 2

)

2

2

p0i−s

1 − f2ik+ −1

1 ∑ 1( s=0

1 ( ∑

1 ∑ 1 s=0

) k+1

(

2

s=0

1 p0i−s + f2ik+ −1

s=0

=

1 ∑ 1

s=0

s=0 1 ∑

ωil−s p0i−s −

1 ∑

ωil−s

s=0

) k+1

p0i−s − f2i−1

) 1 2 1 3 · p0i−s − f2ik+ −1 = O(hk ) · O(hk ) = O(hk ). (

So, if we use the previous relation, we obtain: 1 ωil−1 p0i−1 + ωil p0i − f2ik+ −1 =

1 ∑

1 ωil−s p0i−s − f2ik+ −1

s=0

=( =

1 ∑

1 ωil−s p0i−s − p1i−1 ) + (p1i−1 − f2ik+ −1 )

s=0 O(h3k )

+ O(h2k ) = O(h2k ).

Note that if f has a jump in the interval [ order of the jump, i.e.

xki−1

(11)

□

, ], then the accuracy of the approximation given by xki

p1i−1

is limited by the

1 p1i−1 − f2ik+ −1 = O([f ]),

and proceeding as in (11) we obtain the following error formula: k+1 k l k ωil−1 fi− 1 + ωi fi − f2i−1 = O([f ]).

(12)

We also have Proposition 3. If f is smooth in [xki−2 , xki ] and it has a discontinuity at [xki , xki+1 ] then

ωil−1 =

1+d 2+d

+ O(h) and ωil =

1 2+d

+ O(h) where d = α (f ′ (xki ))2 , l = 0, 1.

(13)

Proof. Using Taylor expansions we obtain the result. □ Please cite this article in press as: F. Aràndiga, Adaptive rational interpolation for point values, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.052.

4

F. Aràndiga / Journal of Computational and Applied Mathematics (

)

–

Proposition 4. If f is smooth in [xki−2 , xki ] and it has a discontinuity at [xki , xki+1 ] then k+1 1 l 0 0 l fˆ2ik+ −1 = ωi−1 pi−1 + ωi pi = f2i−1 + O(hk ),

l = 0, 1.

(14)

Proof. By Proposition 3 we have ωil = O(1) and ωil−1 = O(1), l = 0,1. Then k+1 1 0 l k+1 0 l ˆ k+1 f2ik+ −1 − f2i−1 = ωi−1 (f2i−1 − pi−1 ) + ωi (f2i−1 − pi )

= O(1)O(hk ) + O(1)O(hk ) = O(hk ).

□

3. Rational interpolation of order four k k k k ˆ k+1 We observe that in (10) we use fi− 2 , fi−1 , fi and fi+1 to construct f2i−1 and that

p3i−2 =

1 2

p2i−2 +

1 2

p2i−1 .

where 6 k 3 1 k + fi− + fik , p2i−2 = − fi− 8 2 8 1 8 3 6 1 k p2i−1 = f k + fik − fi+ , 8 i−1 8 8 1 1 k 9 k 9 k 1 k p3i−2 = − fi− + f + f − f . 16 2 16 i−1 16 i 16 i+1

(15)

Then, we have Proposition 5. If f is smooth in [xki−2 , xki+1 ] then 1 k+1 l 2 l 2 4 fˆ2ik+ −1 = ωi−1 pi−2 + ωi pi−1 = f2i−1 + O(hk ),

l = 0, 1.

(16)

Proof. First we prove that 1 ∑

ωil−s p2i−s−1 − p3i−2 =

s=0

=

1 ∑

=

ωil−s p2i−s−1 −

s=0

ωil−s p2i−s−1 −

s=0 1 ∑

1 ∑

1 ∑ 1 s=0

2

s=0

1 p2i−s−1 + f2ik+ −1 1 ∑ 1(

s=0

=

2

s=0

ωil−s −

s=0

1 2

)

1 ∑ 1 s=0

( ) 1 ωil−s p2i−s−1 − f2ik+ −1 −

1 ( ∑

1 ∑ 1

2

2

p2i−s−1

1 − f2ik+ −1

1 p2i−s−1 − f2ik+ −1

1 ∑

ωil−s

s=0

)

) 1 2 3 5 · p2i−s−1 − f2ik+ −1 = O(hk ) · O(hk ) = O(hk ). (

So, if we use the previous relation, we obtain: 1 ωil−1 p2i−2 + ωil p2i−1 − f2ik+ −1 =

1 ∑

1 ωil−s p2i−s−1 − f2ik+ −1

s=0

=(

1 ∑

1 ωil−s p2i−s−1 − p3i−2 ) + (p3i−2 − f2ik+ −1 )

s=0

= O(h5k ) + O(h4k ) = O(h4k ).

□

(17)

But, we also have Proposition 6. Assume that f has a discontinuity at [xki , xki+1 ] and that it is smooth at [xki−2 , xki ]. Then, 1 k+1 l 2 l 2 fˆ2ik+ −1 = ωi−1 pi−2 + ωi pi−1 = f2i−1 + O(1), l = 0, 1.

(18)

Please cite this article in press as: F. Aràndiga, Adaptive rational interpolation for point values, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.052.

F. Aràndiga / Journal of Computational and Applied Mathematics (

)

–

5

Proof. By Proposition 3 we obtain that ωil−s = O(1), s, l = 0,1. Then, k+1 1 2 l k+1 2 l ˆ k+1 f2ik+ −1 − f2i−1 = ωi−1 (f2i−1 − pi−2 ) + ωi (f2i−1 − pi−1 )

= O(1)O(h3k ) + O(1)O(1) = O(1).

□

The results of this section motivate the definition of a new set of weights with improved properties with respect to the accuracy of the resulting interpolant. 4. New weights for rational interpolation In the ENO interpolation [9] we select, between the 2 stencils of 3 points that contains [xki , xki+1 ], the one that minimizes the divided differences of the given data. This selection process leads to a third order interpolant in smooth regions, up to the intervals that contain jump discontinuities, as long as they are ‘‘well separated’’. The idea behind the WENO reconstruction (Liu et al. [2]) is to use the 4 points involved in the ENO reconstruction in order to obtain an approximation of order 4 when f is smooth in [xki−2 , xki+1 ] but also obtain an approximation of order 3 at intervals located next to the one containing a singularity. We achieve this by defining 2,t

2,t

1 2 2 fˆ2ik+ −1 = wi−1 pi−2 + wi pi−1

where

α 2,t ωi2−,ts = ∑2 i−s 2,t , i=0 αi−s

αi2−,ts =

1/2 (ϵ + Ii−s )t

, s = 0, 1,

(19)

and Ii−s is the smoothness indicator of the function f and ϵ is a small positive number introduced to avoid null denominator (in all this paper we take ϵ = 1/h2k as suggested in [10]). Following Jiang and Shu [9] we use Ii−1 = Ii

=

13 12 13 12

k k k 2 (fi− 2 − 2fi−1 + fi ) + k k k 2 (fi− 1 − 2fi + fi+1 ) +

1 4 1 4

k k k 2 (fi− 2 − 4fi−1 + 3fi )

(20)

k k 2 (fi+ 1 − fi−1 )

(21)

as smooth indicators. The weights we propose are inspired on the WENO weights and are specifically designed so that they satisfy the WENO properties. That is, if the function is smooth wi−s = 1/2 + O(h2 ) (Propositions 7 and 8). And if the function has a discontinuity in [xki , xki+1 ] then wi−1 = O(1) and wi = O(h2 ) (Propositions 9 and 10). To obtain this we modify the proposed rational weights 2t of S. Carrato and G. Ramponi ((3) and (4)) [3]. We take α = h− and Ii−s ((20) and (21)) obtaining k

ωi3−,t1 = ωi3,t =

2t t 1 + h− k Ii 2t t t 2 + h− k (Ii + Ii−1 ) 2t t 1 + h− k Ii−1 2t t t 2 + h− h (Ii + Ii−1 )

. .

(22)

Another possibility is to define:

ωi4−,t1 = ωi4,t =

2t t −2t t t 1 + 2h− k Ii + hk Ii Ii−1 2t t t t t 2 + 2h− k (Ii + Ii−1 + Ii Ii−1 ) 2t t −2t t t 1 + 2h− k Ii−1 + hk Ii Ii−1 2t t t t t 2 + 2h− k (Ii + Ii−1 + Ii Ii−1 )

, .

(23)

Now, we have Proposition 7. If f is smooth in [xki−2 , xki+1 ] then ,t ωil− s =

1 2

+ O(h2k ), l = 2, 3, 4.

(24)

Proof. Taking Taylor expansions we obtain Ii−s = (hf ′ (xk2i+−11 ))2 (1 + O(h2k )) and (24). □ We also have Please cite this article in press as: F. Aràndiga, Adaptive rational interpolation for point values, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.052.

6

F. Aràndiga / Journal of Computational and Applied Mathematics (

)

–

Proposition 8. If f is smooth in [xki−2 , xki+1 ] then l,t

l,t

l,t 2 i−1 pi−2

l,t 2 i pi−1

1 k+1 2 0 0 fˆ2ik+ −1 = ωi−1 pi−1 + ωi pi = f2i−1 + O(hk ), l = 2, 3, 4. 1 f2ik+ −1

ˆ

=ω

+ω

1 f2ik+ −1

=

+

O(h4k )

, l = 2, 3, 4.

(25) (26)

Proof. Since we have (24) with the same proofs of Propositions 2 and 5 we obtain (25) and (26). □ With these weights we have Proposition 9. If f is smooth in [xki−2 , xki ] and has a singularity at [xki , xki+1 ] then l,t ,t 2t 2t ωil− 1 = 1 + O(hk ) and ωi = O(hk ) l = 2, 3, 4.

(27)

Proof. Taking Taylor expansions we obtain Ii−1 = (hf ′ (xk2i+−11 ))2 (1 + O(h2k )), Ii = O(1) and (27). □ Proposition 10. Assume that f has a discontinuity at [xki , xki+1 ] and that it is smooth at [xki−2 , xki ]. Then l,1

l,1

l,2

l,2

k+1 1 0 0 1 fˆ2ik+ −1 = ωi−1 pi−1 + ωi pi = f2i−1 + O(hk ), l = 2, 3, 4, k+1 1 3 2 2 fˆ2ik+ −1 = ωi−1 pi−2 + ωi pi−1 = f2i−1 + O(hk ), l = 2, 3, 4. l,t

l,t

2t Proof. We have ωi = O(h2t k ) and ωi−1 = 1 + O(hk ), l = 2, 3, 4. Then l,1

l,1

k+1 k+1 1 0 0 ˆ k+1 f2ik+ −1 − f2i−1 = ωi−1 (f2i−1 − pi−1 ) + ωi (f2i−1 − pi )

=

O(1)O(h1k )

+

O(h2k )O(h1k )

= O(h1k )

and l,2

l,2

k+1 k+1 1 2 2 ˆ k+1 f2ik+ −1 − f2i−1 = ωi−1 (f2i−1 − pi−2 ) + ωi (f2i−1 − pi−1 )

= O(1)O(h3k ) + O(h4k )O(1) = O(h3k ). Observe that in case f has a jump in xd ∈ [ order of the jump, i.e. 1 ∑

xki−1

□

, ] then p2i−s − f2ik −1 = O([f ]), s = 1, 2, and the order of the error is the xki

ωi−s p2i−1−s − f2ik −1 = O([f ]).

(28)

s=0

5. Rational extrapolation of fourth order 1 1 1 1 Let us now define an approximation to f2ik+ −1 with the interpolated values by the polynomials qi−2 (x), qi−1 (x) and qi (x) (see k (5)). The values of these polynomials at x2i−1 are given by the following expressions:

1 k 3 k p1i−2 = − fi− + fi− , 2 2 2 1 1 k 1 p1i−1 = fi−1 + fik , 2 2 3 k 1 k 1 pi = f − fi+1 . 2 i 2

(29)

If we consider the previous expressions and p3i−2 (15), then we have: p3i−2 =

1 8

p1i−2 +

6 8

p1i−1 +

1 8

p1i .

(30)

In Sections 3 and 4 we have used a weighted average of the evaluation at xk2i+−11 ∈ [xki−1 , xki ] of the two polynomials that interpolate the points {xki−2 , xki−1 , xki } and {xki−1 , xki , xki+1 }. In this section we use three polynomials. In two of them, the polynomial is evaluated at a point xk2i+−11 outside the interpolation nodes {xki−2 , xki−1 } and {xki , xki+1 }. For this reason, and to distinguish it from the previous case, we call it extrapolation. 1 Now, we define an approximation to f2ik+ −1 as follows. 5,t

5,t

5 ,t

1 1 1 1 fˆ2ik+ −1 = ωi−2 pi−2 + ωi−1 pi−1 + ωi pi ,

(31)

Please cite this article in press as: F. Aràndiga, Adaptive rational interpolation for point values, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.052.

F. Aràndiga / Journal of Computational and Applied Mathematics (

)

–

7

where

ωi5−,t2 = ωi5−,t1 = ωi5,t =

2t 2 −4t 2 2 2 1 + h− k (Ji−1 + Ji ) + hk Ji−1 Ji

−2t

8 + 2hk

(Ji2−2

4t 2 2 2 2 2 2 + Ji2−1 + Ji2 ) + h− k (Ji−1 Ji + Ji−2 Ji + Ji−2 Ji−1 )

2t 2 −4t 2 2 2 6 + 6h− k (Ji−2 + Ji ) + 6hk Ji−2 Ji

−4t 2 2 2t 2 2 2 2 2 2 2 8 + 2h− k (Ji−2 + Ji−1 + Ji ) + hk (Ji−1 Ji + Ji−2 Ji + Ji−2 Ji−1 ) 2t 2 −4t 2 2 2 1 + h− k (Ji−2 + Ji−1 ) + hk Ji−2 Ji−1 2t 2 −4t 2 2 2 2 2 2 2 2 8 + 2h− k (Ji−2 + Ji−1 + Ji ) + hk (Ji−1 Ji + Ji−2 Ji + Ji−2 Ji−1 )

, ,

(32)

,

k k 2 and Ji−s = (fi− s+1 − fi−s ) , s = 0, 1, 2, t = 1, 2. Now, we have

Proposition 11. If f is smooth in [xki−2 , xki+1 ], then 5,t

5,t

5 ,t

k+1 1 1 4 1 1 fˆ2ik+ −1 = ωi−2 pi−2 + ωi−1 pi−1 + ωi pi = f2i−1 + O(hk ).

(33)

Proof. First, we assume that t = 1. If f (x) is smooth in [xki−2 , xki+1 ] taking Taylor expansions we can verify the following relations:

ωi5−,12 − C−2 = O(h2k ), where C−2 = C0 = 2 ∑

1 8

6 . 8

and C−1 =

ωi5−,1s p1i−s − p3i−2 =

s=0

2 ∑

ωi5−,1s p1i−s −

2 ∑

C−s p1i−s

s=0

2 ∑

1 ωi5−,1s p1i−s − f2ik+ −1

s=0

2 ∑

2 ∑

C−s −

s=0 2 ∑

2 ∑

C−s p1i−s

s=0

1 ωi5−,1s (p1i−s − f2ik+ −1 ) −

=

2 ∑

1 C−s (p1i−s − f2ik+ −1 )

s=0

s=0 2 ∑

ωi5−,1s

s=0

1 +f2ik+ −1

=

(34)

Then,

s=0

=

ωi5,1 − C0 = O(h2k ),

ωi5−,11 − C−1 = O(h2k ),

5,1

1 2 2 (ωi−s − C−s )(p1i−s − f2ik+ −1 ) = O(hk )O(hk )

s=0

= O(h4k ).

(35)

Now, we assume that t = 2. If f (x) is smooth in [xki−2 , xki+1 ] taking Taylor expansions we can verify the following relations:

ωi5−,22 =

1 8

+h

f2i′′−1 (f2i′ −1 )3 ′

2(f2i−1 )4 + 2

+ O(h2k ), ωi5−,21 =

6 8

+ O(h2k ), ωi5,2 =

1 8

−h

f2i′′−1 (f2i′ −1 )3 2(f2i′ −1 )4 + 2

+ O(h2k )

(36)

where f2i′ −1 = f ′ (xk2i+−11 ) and f2i′′−1 = f ′′ (xk2i+−11 ). We also have 3 ′′ 1 2 p1i−s − f2ik+ −1 = − f2i−1 hk , s = 0, 2. 8

(37)

Then 2 ∑

ωi5−,2s p1i−s − p3i−2 =

s=0

2 ∑

5,t

1 2 2 (ωi−s − C−s )(p1i−s − f2ik+ −1 ) = O(hk )O(hk )

s=0

+hk =

f2i′′−1 (f2i′ −1 )3 2 3 ′′ 3 h2k f2i′′−1 − hk h f +2 8 2(f2i′ −1 )4 + 2 k 8 2i−1

f2i′′−1 (f2i′ −1 )3 ′

2(f2i−1

O(h4k )

.

)4

(38)

Please cite this article in press as: F. Aràndiga, Adaptive rational interpolation for point values, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.052.

8

F. Aràndiga / Journal of Computational and Applied Mathematics (

)

–

Table 1 2,1 4,1 1 k l l 1 ˆ k+1 Errors elk = max{|f2ik+ −1 − f2i−1 |, xi ∈ [0, 0.25]} and orders of approximation log2 (ek /ek−1 ) for the weights ωi−s , (3), ωi−s , (19) and ωi−s , (23). 2,1

1 0 1 1 0 fˆ2ik+ −1 = ωi−1 pi−1 + ωi pi

2,1

4,1

1 0 0 fˆ2ik+ −1 = ωi−1 pi−1 + ωi pi

/

e1k

log2 (e1k

2.3581e−05 5.9186e−06 1.4825e−06 3.7100e−07 9.2795e−08 2.3204e−08 5.8018e−09 1.4505e−09

1.9943e+00 1.9972e+00 1.9986e+00 1.9993e+00 1.9996e+00 1.9998e+00 1.9999e+00 2.0000e+00

e1k−1 )

4,1

1 0 0 fˆ2ik+ −1 = ωi−1 pi−1 + ωi pi

/

e2k

log2 (e2k

2.3297e−05 5.8833e−06 1.4781e−06 3.7045e−07 9.2726e−08 2.3196e−08 5.8007e−09 1.4504e−09

1.9855e+00 1.9928e+00 1.9964e+00 1.9982e+00 1.9991e+00 1.9996e+00 1.9998e+00 1.9999e+00

e2k−1 )

e4k

log2 (e4k /e4k−1 )

2.3168e−05 5.8673e−06 1.4762e−06 3.7020e−07 9.2695e−08 2.3192e−08 5.8003e−09 1.4503e−09

1.9814e+00 1.9908e+00 1.9954e+00 1.9977e+00 1.9989e+00 1.9994e+00 1.9997e+00 1.9999e+00

If we use (35) or (38), we obtain the order of the technique: 2 ∑

1 ωi5−,ts p1i−s − f2ik+ −1 = (

s=0

2 ∑

1 ωi5−,ts p1i−s − p3i−2 ) + (p3i−2 − f2ik+ −1 )

s=0

= O(h4k ) + O(h4k ) = O(h4k ).

(39)

□

Proposition 12. Assume that f has a discontinuity at [xki , xki+1 ] and that it is smooth at [xki−2 , xki ]. Then 5 ,t

5,t

5,t

k+1 1 1 1 2 1 fˆ2ik+ −1 = ωi−2 pi−2 + ωi−1 pi−1 + ωi pi = f2i−1 + O(hk ), t = 1, 2.

(40) 5 ,t

5 ,t

5 ,t

Proof. Using Taylor expansions we obtain that ωi−2 = 1/7 + O(h2k ), ωi−1 = 6/7 + O(h2k ) and ωi 5,t

5,t

= O(h2t k ). Then

5 ,t

k+1 1 k+1 k+1 1 1 1 ˆ k+1 f2ik+ −1 − f2i−1 = ωi−2 (f2i−1 − pi−2 ) + ωi−1 (f2i−1 − pi−1 ) + ωi (f2i−1 − pi )

= O(1)O(h2k ) + O(1)O(h2k ) + O(h2k )O(1) = O(h2k ).

□ 5,t

5,t

Proposition 13. if f has a discontinuity at [xki−1 , xki ] and is smooth elsewhere then ωi−2 = 1/2 + O(h2k ), ωi−1 = O(h2t k ) and ωi5,t = 1/2 + O(h2k ). Proof. Using Taylor expansions. □

6. Numerical results 6.1. Order of reconstructions In this experiment we validate the theoretical results obtained in the previous sections on the order of approximation of our techniques. We consider the function

{ f (x) =

ex x 1 + ex N

if if

0 ≤ x ≤ 0.5 0.5 < x ≤ 1.

(41) N

1 From f k = {fik }i=k0 , Nk = 2k , k = 6 : 13, we obtain an approximation fˆ k+1 to f k+1 = {fik+1 }i=k+ 0 using the different strategies presented in this paper. First, we want to see if we obtain the expected order when f is smooth. 1 k l l ˆ k+1 We evaluate ek = max{|f2ik+ −1 − f2i−1 |, xi ∈ [0, 0.25]} (note that f is smooth at [0, 0.25]) and log2 (ek /ek−1 ) with ∑ l n l l ˆf k+1 = 1 ωl pr 2i−1 s=0 i−s i−s−r /2 . Note that if ek ≈ hk then log2 (ek /ek−1 ) ≈ n. We compare the results that we obtain when we use the rational weights of S. Carrato and G. Ramponi [3], ωi1−s (3), the 2,t 4,t WENO weights, ωi−s (19), and those presented in this paper, ωi−s (23). We perform the experiment for r = 0 (Table 1) and for r = 2 (Table 2). We observe that we obtain the expected orders of approximation. That is, order 2 for r = 0 (Propositions 2 and 8) and 4 for r = 2 (Propositions 5 and 8). Now we evaluate the error in [0, 0.5] where the singularity is precisely the end-point of the considered region. Note that when we use linear reconstruction techniques with 4 points the accuracy of the approximation at the interval [0.5 − hk , 0.5] is affected by the discontinuity and the approximation order is O(1). We perform the experiment for r = 0 (Table 3) and for r = 2 (Table 4).

Please cite this article in press as: F. Aràndiga, Adaptive rational interpolation for point values, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.052.

F. Aràndiga / Journal of Computational and Applied Mathematics (

)

–

9

Table 2 2 ,2 4,2 1 k l l 1 ˆ k+1 Errors elk = max{|f2ik+ −1 − f2i−1 |, xi ∈ [0, 0.25]} and orders of approximation log2 (ek /ek−1 ) for the weights ωi−s , (3), ωi−s , (19) and ωi−s , (23). 2,2

1 1 2 1 2 fˆ2ik+ −1 = ωi−1 pi−2 + ωi pi−1

2,2

4,2

1 2 2 fˆ2ik+ −1 = ωi−1 pi−2 + ωi pi−1

/

e1k

log2 (e1k

1.4412e−09 8.9163e−11 5.5445e−12 3.4565e−13 2.1576e−14 1.3476e−15 8.4200e−17 5.2616e−18

4.0147e+00 4.0074e+00 4.0037e+00 4.0018e+00 4.0009e+00 4.0005e+00 4.0003e+00 4.0001e+00

e1k−1 )

4,2

1 2 2 fˆ2ik+ −1 = ωi−1 pi−2 + ωi pi−1

/

e2k

log2 (e2k

1.4412e−09 8.9163e−11 5.5445e−12 3.4565e−13 2.1576e−14 1.3476e−15 8.4200e−17 5.2616e−18

4.0147e+00 4.0074e+00 4.0037e+00 4.0018e+00 4.0009e+00 4.0005e+00 4.0003e+00 4.0001e+00

e2k−1 )

e4k

log2 (e4k /e4k−1 )

1.4412e−09 8.9163e−11 5.5445e−12 3.4565e−13 2.1576e−14 1.3476e−15 8.4200e−17 5.2616e−18

4.0147e+00 4.0074e+00 4.0037e+00 4.0018e+00 4.0009e+00 4.0005e+00 4.0003e+00 4.0001e+00

Table 3 2,1 4,1 1 k l l 1 ˆ k+1 Errors elk = max{|f2ik+ −1 − f2i−1 |, xi ∈ [0, 0.5]} and orders of approximation log2 (ek /ek−1 ) for the weights ωi−s , (3), ωi−s , (19) and ωi−s , (23). 2,1

1 0 1 1 0 fˆ2ik+ −1 = ωi−1 pi−1 + ωi pi

2,1

4,1

1 0 0 fˆ2ik+ −1 = ωi−1 pi−1 + ωi pi

/

e1k

log2 (e1k

e1k−1 )

4.7649e−05 3.1222e−05 1.7471e−05 9.2018e−06 4.7176e−06 2.3880e−06 1.2013e−06 6.0248e−07

6.0985e−01 8.3762e−01 9.2501e−01 9.6384e−01 9.8222e−01 9.9117e−01 9.9559e−01 9.9783e−01

4,1

1 0 0 fˆ2ik+ −1 = ωi−1 pi−1 + ωi pi

/

e2k

log2 (e2k

e2k−1 )

7.7157e−03 3.8827e−03 1.9473e−03 9.7512e−04 4.8792e−04 2.4405e−04 1.2205e−04 6.1030e−05

9.9074e−01 9.9559e−01 9.9783e−01 9.9892e−01 9.9949e−01 9.9971e−01 9.9986e−01 9.9993e−01

e4k

log2 (e4k /e4k−1 )

7.7153e−03 3.8827e−03 1.9473e−03 9.7512e−04 4.8792e−04 2.4405e−04 1.2205e−04 6.1030e−05

9.9066e−01 9.9559e−01 9.9783e−01 9.9892e−01 9.9949e−01 9.9971e−01 9.9986e−01 9.9993e−01

Table 4 2,2 4,2 1 k l l 1 ˆ k+1 Errors elk = max{|f2ik+ −1 − f2i−1 |, xi ∈ [0, 0.5]} and orders of approximation log2 (ek /ek−1 ) for the weights ωi−s , (3), ωi−s , (19) and ωi−s , (23). 2,2

2,4

2,2

/

e1k

log2 (e1k

e1k−1 )

6.0912e−02 6.1395e−02 6.1638e−02 6.1759e−02 6.1820e−02 6.1851e−02 6.1866e−02 6.1874e−02

−1.1384e−02 −5.6954e−03 −2.8449e−03 −1.4290e−03 −7.0709e−04 −3.6072e−04 −1.7313e−04 −8.6564e−05

2,4

1 2 2 fˆ2ik+ −1 = ωi−1 pi−2 + ωi pi−1

1 2 2 fˆ2ik+ −1 = ωi−1 pi−2 + ωi pi−1

1 1 2 1 2 fˆ2ik+ −1 = ωi−1 pi−2 + ωi pi−1

/

e2k

log2 (e2k

e2k−1 )

2.1805e−07 2.8540e−08 3.6467e−09 4.6076e−10 5.7902e−11 7.2569e−12 9.0830e−13 1.1361e−13

2.9336e+00 2.9683e+00 2.9845e+00 2.9923e+00 2.9962e+00 2.9981e+00 2.9990e+00 2.9995e+00

e4k

log2 (e4k /e4k−1 )

−4.7491e+00 2.8971e−08 3.6734e−09 4.6242e−10 5.8005e−11 7.2633e−12 9.0870e−13 1.1364e−13

2.9578e+00 2.9794e+00 2.9898e+00 2.9949e+00 2.9975e+00 2.9987e+00 2.9994e+00 2.9997e+00

Table 5 1 k ˆ k+1 Errors elk = max{|f2ik+ −1 − f2i−1 |, xi ∈ [0, a]} and orders of approximation 5 ,2

5,2

5 ,2

1 2 2 2 log2 (elk /elk−1 ) for the approximations fˆ2ik+ −1 = ωi−2 pi−2 + ωi−1 pi−1 + ωi pi

with weights ω

5 ,2 i−s

(32).

Error in ∈ [0, 0.25]

Error in [0, 0.5]

e5k

log2 (e5k /e5k−1 )

e5k

log2 (e5k /e5k−1 )

6.2523e−09 3.9418e−10 2.4741e−11 1.5495e−12 9.6946e−14 6.0623e−15 3.7899e−16 2.3690e−17

3.9875e+00 3.9939e+00 3.9970e+00 3.9985e+00 3.9993e+00 3.9996e+00 3.9998e+00 3.9999e+00

1.2653e−05 3.2160e−06 8.1070e−07 2.0351e−07 5.0984e−08 1.2759e−08 3.1915e−09 7.9807e−10

1.9761e+00 1.9880e+00 1.9940e+00 1.9970e+00 1.9985e+00 1.9992e+00 1.9996e+00 1.9998e+00

We obtain the expected orders: 1 for r = 0 (Propositions 4 and 10), 0 for r = 2 and the weights of S. Carrato and G. l,2 Ramponi [3] ωi1 (Proposition 6) and 3 for r = 2 and the new weights ωi , l = 2,4 (Proposition 10). Now, we consider the rational extrapolation. We observe in Table 5 that we obtain the expected order 4 (Proposition 11) when f is smooth (interval [0, 0.25]) and 2 (Proposition 12) when we are close to a discontinuity (interval [0, 0.5]). Please cite this article in press as: F. Aràndiga, Adaptive rational interpolation for point values, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.052.

10

F. Aràndiga / Journal of Computational and Applied Mathematics (

)

–

1 r r Fig. 1. Reconstruction with fˆ2ik+ −1 = ωi−1 pi−1−r /2 + ωi pi−r /2 . Left, r = 0. Right, r = 2. From top to bottom: ωi−1 = ωi =

1 , 2

2,r /2+1

ωi1−s , ωi−s

4,r /2+1

and ωi−s

.

Please cite this article in press as: F. Aràndiga, Adaptive rational interpolation for point values, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.052.

F. Aràndiga / Journal of Computational and Applied Mathematics (

)

–

11

5 ,t 1 5 ,t 1 5,t 1 1 Fig. 2. fˆ2ik+ −1 = ωi−2 pi−2 + ωi−1 pi−1 + ωi pi . Left t = 1. Right t = 2.

i , Fig. 3. Function (42). Left {f ( 10

j ) 10 . 10 i,j=0

}

i Right: {f ( 40 ,

j ) 40 . 40 i,j=0

}

6.2. Limit function in 1D

We take the same function (41). We apply recursively the different interpolatory techniques from the resolution level 3 7 k = 3 ({f (x3i )}2i=0 ) to obtain, an approximation to the data at the resolution level k = 7 ({f (x7i )}2i=0 ). The results we obtain are shown in Figs. 1 and 2. The line convention in the figures is as follows.

• Circles: function in resolution level k = 3. • Points: function in resolution level k = 7. • Dashed Line: reconstruction. 1 r r We take the reconstruction fˆ2ik+ −1 = ωi−1 pi−1−r /2 + ωi pi−r /2 (Fig. 1). When ωi−1 = ωi =

1 2

we obtain an approximation of

order 2 for r = 0 and O(1) when r = 2. 2,r /2+1 4,r /2+1 If we take r = 0 and the nonlinear weights (ωi1−s , ωi−s , ωi−s ) we always obtain an approximation of order 1 but k k the different behavior of the weights near the discontinuity (Propositions 3 and 9) causes that the influence of fi− 1 and fi is different, obtaining what we see in the figure. If r = 2 we obtain an approximation of order 1 when we take the weights of S. Carrato and G. Ramponi [3] and 3 with the weights of Section 4. 5 ,t 1 5,t 1 5,t 1 1 In Fig. 2 we show the result we obtain when we take fˆ2ik+ −1 = ωi−2 pi−2 + ωi−1 pi−1 + ωi pi . The contribution of the polynomial of the interval that contains the discontinuity is smaller when t = 2 (Proposition 13) obtaining a sharper reconstruction in this case. Please cite this article in press as: F. Aràndiga, Adaptive rational interpolation for point values, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.052.

12

F. Aràndiga / Journal of Computational and Applied Mathematics (

)

–

2,r /2+1

Fig. 4. Reconstruction with wi−1 pri−1−r /2 + wi pri−r /2 . Left r = 0. Right r = 2. From Top to bottom wi1−s , wi−s

4,r /2+1

wi−s

.

6.3. Reconstruction in 2D In this case f (x, y) is defined as (42)

⎧ sin(π x) sin(π y) ⎪ ⎨ − sin(π x) sin(π y) + 2 f (x, y) = ⎪− sin(π x) sin(π y) + 2 ⎩ − sin(π x) sin(π y) + 4

if if if if

0 ≤ x ≤ .5; 0 ≤ y ≤ .5 .5 < x ≤ 1; 0 ≤ y ≤ .5 0 ≤ x ≤ .5; .5 < y ≤ 1 .5 < x ≤ 1; .5 < y ≤ 1,

(42)

i i 40 and from {f ( 10 , 10 )}10 i,j=0 we obtain an approximation to {f ( 40 , 40 )}i,j=0 (Fig. 3) using the different reconstructions of this paper via tensor product. That is, first the 1D technique is applied by rows, then by columns and so on. 1 r r In Fig. 4 we show what we obtain with the reconstructions fˆ2ik+ −1 = ωi−1 pi−1−r /2 + ωi pi−r /2 . When we take r = 0 and j

j

2,r /2+1

4,r /2+1

the nonlinear weights ( ωi1−s , ωi−s and ωi−s ) we always obtain an approximation of order 1 close to the jump. The k k different behavior of the weights near the discontinuity (Propositions 3 and 9) causes that the influence of fi− 1 and fi is different obtaining what we see in the figure.

Please cite this article in press as: F. Aràndiga, Adaptive rational interpolation for point values, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.052.

F. Aràndiga / Journal of Computational and Applied Mathematics (

5 ,t

5,t

)

–

13

5,t

Fig. 5. Reconstruction with wi−2 p1i−2 + wi−1 p1i−1 + wi p1i . Left t = 1. Right t = 2.

If r = 2 we obtain near the jump an approximation of order 1 when we take the weights of S. Carrato and G. Ramponi [3] (we can observe the Gibbs phenomenon) and 3 with the weights of Section 4. 5,t 1 5,t 1 5,t 1 1 In Fig. 5 we show the result we obtain when we take fˆ2ik+ −1 = ωi−2 pi−2 + ωi−1 pi−1 + ωi pi . The contribution of the polynomial of the interval that contains the discontinuity is smaller when t = 2 (Proposition 13) obtaining a sharper reconstruction in this case. References [1] A. Harten, B. Engquist, S. Osher, S. Chakravarthy, Uniformly high order accurate essentially non-oscillatory schemes III, J. Comput. Phys. 71 (1987) 231–303. [2] X.-D. Liu, S. Osher, T. Chan, Weighted Essentially Non-oscillatory Schemes, J. Comput. Phys. 115 (1994) 200–212. [3] S. Carrato, G. Ramponi, Interpolation of the DC component of coded images using a rational filter, in: Proc. Fourth IEEE Intern. Conf. on Image Processing - 97, pp. 26–29, S. Barbara, CA, Oct. 1997. [4] R. Castagno, S. Marsi, G. Ramponi, A simple algorithm for the reduction of blocking artifacts in images and its implementation, IEEE Trans. Consum. Electron. 44 (3) (1998) 1062–1070. [5] R. Castagno, G. Ramponi, A rational filter for the removal of blocking artifacts in image sequences coded at low bitrate, in: Proc. Eighth European Signal Processing Conf. EUSIPCO-96, Trieste, Italy, Sept. 10-13 1996. [6] G. Ramponi, Image Processing using Rational Functions, in: Proc. Second IEEE Intern. Conf. on Image Processing - 95, pp. 22–25, Washington, DC, USA, Oct 1995. [7] A. Harten, Discrete multiresolution analysis and generalized wavelets, J. Appl. Num. Math. 12 (1993) 153–193. [8] A. Harten, Multiresolution representation of data: a general framework, SIAM J. Numer. Anal. 33 (3) (1996) 1205–1256. [9] G.-S. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1) (1996) 202–228. [10] F. Aràndiga, A.M. Belda, P. Mulet, Point-value WENO multiresolution applications to stable image compression, J. Sci. Comput. 43 (2) (2010) 158–182.

Please cite this article in press as: F. Aràndiga, Adaptive rational interpolation for point values, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.08.052.