I want to generate a small random perturbation on a flat surface defined on a square domain, e.g. 1 + 0.05 rand(x,y), which can be used as initial condition with periodic boundary conditions in `NDSolve`. In other word, rand(x,y) is a special pseudo-random function distributed in the interval (-1,1) which has periodicity at the sides of the square-domain [0,L]*[0,L], that is, rand(0,y)=rand(L,y) and rand(x,0)=rand(x,L). In fact, this question is very similar to [this question][1]. 

**However, after I went through this post in details, all comment, voted answers and even "accepted" answer as well as the answer given by the author, I found this question is still unsolved.**

 
Based on [acl's comment][2] and [Vitaliy's answer][3], I have the following working example with two different initial conditions.

    Off[NDSolve::mxsst];
    L = 4*\[Pi];
    c = -(1/20);
    tmax = 1000;
    ini1 = Interpolation@Flatten[Table[{{x, y}, 1 + c*RandomReal[{-1, 1}]}, {x, 0, L + 1}, {y, 0, L + 1}], 1]; (*initial condition given by acl*)
    cond3D = First[h /. NDSolve[{D[h[x, y, t], t] + 
    Div[h[x, y, t]^3*Grad[Laplacian[h[x, y, t], {x, y}], {x, y}], {x, y}] + 
    Div[h[x, y, t]^3*Grad[h[x, y, t], {x, y}], {x, y}] == 0,
    h[0, y, t] == h[L, y, t],
    h[x, 0, t] == h[x, L, t],
    h[x, y, 0] == ini1[x, y]
     },
    h, {x, 0, L}, {y, 0, L}, {t, 0, tmax},
    Method -> {"BDF", "MaxDifferenceOrder" -> 1}, 
    MaxStepFraction -> 1/50]]

I have the following error:
>NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent.


If I use the other initial condition given by [@Vitaliy][4] just with a slight modification:

    ini2 = 1 + c*BSplineFunction[RandomReal[1, {30, 30, 1}], SplineClosed -> True][##] &;

where, because I want a small perturbation so a small parameter `c`. Just replace the initial condition as

    h[x, y, 0] == ini2[x, y]
 The `NDSlove` then eats up the memory and spit a lot a error related to data.

I realized that [Vitaliy's answer][5] is wonderful in `SplineClosed -> True` that makes I can use it with periodic boundary conditions in `NDSolve`. But another confused issue is when I try, `ini2[9, 0] === ini2[9, L] // FullSimplify` using my `ini2`, MMA gives **False**.I know @Vitaliy produced a disturbed surface in the domain[0, 1]* [0, 1], but how to get such a surface in a square domain, [0, L]*[0, L].

How can I generate a random disturbance consisting with periodic boundary condition meanwhile. If anyone can give such an initial condition that has something like the `SplineClosed->True` feature that can be fed in `NDSolve`, that would be truly awesome! I admit I am being very silly when I try to use RandomReal[] in this way.


  [1]: http://mathematica.stackexchange.com/questions/15726/small-random-disturbance-of-a-flat-surface
  [2]: http://mathematica.stackexchange.com/questions/15726/small-random-disturbance-of-a-flat-surface#comment46466_15726
  [3]: http://mathematica.stackexchange.com/questions/15726/small-random-disturbance-of-a-flat-surface#comment46466_15726
  [4]: http://mathematica.stackexchange.com/questions/15726/small-random-disturbance-of-a-flat-surface/15728#15728
  [5]: http://mathematica.stackexchange.com/questions/15726/small-random-disturbance-of-a-flat-surface/15728#15728