Re: [ferret_users] boundary condition for smoothing (@SPZ and friends)

["William S. Kessler" <William.S.Kessler@xxxxxxxx>]

Great suggestion! Carrying the filtering to the edges of the grid would 
not be too hard to write a script for. (It would be a little harder to 
do exactly on a non-uniform grid, such as we often have in the 
meridional direction in models). Filtering across interior missing 
points is less obvious: suppose there is a single missing gridpoint. 
Then the algorithm discussed here would lead to assigning different 
values to that gridpoint to implement the filter from either side. In 
theory that is not a problem but in practice I can't think of a way to 
do that using Ferret commands. It would have to be done in fortran as 
an external function, I think. But I hope someone figures it out as 
this is a much-needed addition.

Billy K

On Sep 2, 2005, at 11:54 AM, Ryo Furue wrote:

> Hello Ferret users,
>
> I noticed that smoothing transformations (@SPZ and friends) leave
> endpoints undefined. For example,
>
>   yes? let pi = 4 * atan(1)
>   yes? let func sin(pi*x[x=-1:1:0.2])
>   yes? plot/line/symbol/hlimits=-1.2:1.2    func
>   yes? plot/line/symbol/hlimits=-1.2:1.2/ov func[x=@spz]
>
> Even though func has datapoints at the edges (x = -1 and x = 1),
> func[x=@spz] has missing values there.  I'm wondering if
> those transformations can be "improved", at least in simple cases.
>
> Let me take @SPZ as an example.  @SPZ is a 1-2-1 weighted moving
> average:
>
>     g(i) = [f(i-1) + 2 f(i) + f(i+1)] / 4
>
> where f is defined for i = 1, 2, . . ., N.  What to do with g(1)
> and g(N), for which f(0) and f(N+1) would be required?
>
> One solution is to define g(1) == f(1) and g(N) == f(N).
> But, I think a better solution is to make the smoothing behave
> like diffusion with a no-flux boundary condition(*).  That is,
> we define extra gridpoints at i = 0 and i = N + 1 and
> define the values of f(i) there as f(0) == f(1) and f(N+1) == f(N).
> And then we compute g(1) and g(N) using the formula above:
>
>     g(1) = [f(0) + 2 f(1) + f(2)] / 4
>          = [f(1) + 2 f(1) + f(2)] / 4
>          = [       3 f(1) + f(2)] / 4
>
> and likewise for g(N).
>
> We could apply the same no-flux condition around "interior"
> missing datapoints, too.
>
> Advantages of this approach are that averages are conserved:
> g[i=1:N@sum] == f[i=1:N@sum]  (This property corresponds to
> the conservation of heat in the diffusion equation with no-flux
> boundary condition.), and that it matches our intuitive
> understanding of "mixing" neighboring values.
>
> I haven't given a serious thought to other lowpass filters; the above
> consideration may or may not apply to them.
>
> A practical reason for proposing this is that I don't like losing
> vectors or shading near the boundaries when smoothing is applied.
>
> Regards,
> Ryo
> ======================================
> (*)In fact, the 1-2-1 moving average can be cast into the diffusion
> equation:
>
> g(i) = [f(i-1) + 2 f(i) + f(i+1)] / 4
> <-->
> g(i) - f(i) = [f(i-1) - 2 f(i) + f(i+1)] / 4
> <-->
> [g(i) - f(i)] / delta_t = kappa [f(i-1) - 2 f(i) + f(i+1)] / 
> (delta_x)^2
>
> with a suitable choice of delta_t, kappa, and delta_x.
> This is a finite-difference form of the diffusion equation
> with g(i) being the value of the next timestep.
>
> I think the 1-4-6-4-1 smoothing can be likewise cast into
> a diffusion equation with biharmonic diffusion.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
William S. Kessler
NOAA / Pacific Marine Environmental Laboratory
7600 Sand Point Way NE
Seattle WA 98115 USA

william.s.kessler@noaa.gov
Tel: 206-526-6221
Fax: 206-526-6744
Home page: http://www.pmel.noaa.gov/~kessler