Small random disturbance of a flat surface

Question

I am trying to create an initial condition which is:

1 + 0.05 rand(x,y) Here rand is a pseudorandom function distributed in the interval (-1,1). This surface represents a random disturbance that I would like to use as an initial condition for PDEs in NDSolve.

I assume I am being very silly when I try to use RandomReal[] as my random number generator for my random disturbance. How should I proceed with this.

L = 100;
Plot3D[
       1 - 0.05 (Cos[2 π x/L] + Sin[2 π x/L]) Cos[2 π y/L] RandomReal[],
       {x, 0, L}, {y, 0, L}
       ]

Obviously, this is wrong as this still retains the underlying Cos/Sin curve. How should I go about creating a random disturbance? $\delta\varepsilon\pi$

Working example:

$HistoryLength = 0;
Needs["VectorAnalysis`"]
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Clear[Eq0, EvapThickFilm, h, Bo, ε, K1, \[Delta], Bi, m, r]
Eq0[h_, {Bo_, ε_, K1_, δ_, Bi_, m_, r_}] := D[h, t] + 
    Div[-h^3 Bo Grad[h] + 
      h^3 Grad[Laplacian[h]] + (δ h^3)/(Bi h + K1)^3 Grad[h] + 
      m (h/(K1 + Bi h))^2 Grad[h]] + ε/(
    Bi h + K1) + (r) D[D[(h^2/(K1 + Bi h)), x] h^3, x] == 0;
SetCoordinates[Cartesian[x, y, z]];
EvapThickFilm[Bo_, ε_, K1_, δ_, Bi_, m_, r_] := 
  Eq0[h[x, y, t], {Bo, ε, K1, δ, Bi, m, r}];
TraditionalForm[
  EvapThickFilm[Bo, ε, K1, δ, Bi, m, r]];


L = 2*92.389; TMax = 3100*100;
Off[NDSolve::mxsst];
Clear[Kvar];
Kvar[t_] :=  Piecewise[{{1, t <= 1}, {2, t > 1}}]
(* Ktemp = Array[0.001 + 0.001 #^2 &, 13] *)
hSol = h /. NDSolve[{
     (*Bo,ε,K1,δ,Bi,m,r*)

     EvapThickFilm[0.003, 0, 1, 0, 1, 0.025, 0],
     h[0, y, t] == h[L, y, t],
     h[x, 0, t] == h[x, L, t],
     (*h[x,y,0] == 1.1+Cos[x] Sin[2y] *)

     h[x, y, 0] == BSplineFunction[RandomReal[1, {30, 30, 1}]]
     },
    h,
    {x, 0, L},
    {y, 0, L},
    {t, 0, TMax},
    Method -> {"BDF", "MaxDifferenceOrder" -> 1},
    MaxStepFraction -> 1/50
    ][[1]]

With the B-spline as suggested by Vitaliy Kaurov in the answer below, I have the following error:

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>

ReplaceAll::reps: {(h^(0,0,1))[x,y,t]-0.009 h[x,y,t]^2 (h^(0,1,0))[x,y,t]^2-(0.05 h[x,y,t]^2 (h^(0,1,0))[x,y,t]^2)/(1+h[x,y,t])^3+<<13>>+h[x,y,t]^3 ((h^(0,4,0))[x,y,t]+(h^(2,2,0))[x,y,t])+3 h[x,y,t]^2 (h^(1,0,0))[x,y,t] ((h^(1,2,0))[x,y,t]+(h^(3,0,0))[x,y,t])+h[x,y,t]^3 ((h^(2,2,0))[x,y,t]+(h^(4,0,0))[x,y,t])==0,h[0,y,t]==h[184.778,y,t],h[<<1>>]==<<1>>,h[x,y,0]==BSplineFunction[{{0.,1.},{0.,1.}},<>]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

Answer 1

I would use splines - it is very easy:

f = BSplineFunction[RandomReal[1, {30, 30, 1}], SplineClosed -> True]
Plot3D[f[x, y], {x, 0, 1}, {y, 0, 1}, 
       ColorFunction -> "DarkRainbow", Mesh -> All, MeshStyle -> Opacity[.2]]

SplineClosed -> True makes sure you can use it with periodic boundary conditions in NDSolve. This is to show that indeed surface has periodic boundary conditions:

f[.7, 0] == f[.7, 1]

True

Manipulate[
 Plot[{f[x, 0], f[x, y]}, {x, 0, 1}, PlotRange -> {0, 1}, Filling -> {1 -> {2}}],
 {y, 0, 1}]

Answer 2

For random disturbances that retain some smoothness, I turn to Perlin noise:

dot2 = With[{grad = Most[Tuples[{1, -1, 0}, {2}]]},
            Compile[{{gradIdx, _Integer}, {x, _Real}, {y, _Real}},
                    {x, y}.grad[[gradIdx + 1]]]];

fade = Compile[{{t, _Real}}, t*t*t/(3.*t*(t - 1.) + 1.), RuntimeAttributes -> {Listable}];

lerp = Compile[{{x, _Real}, {y, _Real}, {t, _Real}}, (1. - t)*x + t*y];

perlin2D = With[{perms = Apply[Join, ConstantArray[RandomSample[Range[0, 15]], 2]]},
                Compile[{{x, _Real}, {y, _Real}},
                        Module[{xi, yi, xa, ya, u, v, g00, g10, g01, g11},

                               xi = Floor[x]; yi = Floor[y];
                               xa = x - xi; ya = y - yi;
                               xi = Mod[xi, 16] + 1; yi = Mod[yi, 16] + 1;

                               u = fade[xa]; v = fade[ya];

                               g00 = Mod[perms[[perms[[xi]] + yi]], 8];
                               g10 = Mod[perms[[perms[[xi + 1]] + yi]], 8];
                               g01 = Mod[perms[[perms[[xi]] + yi + 1]], 8];
                               g11 = Mod[perms[[perms[[xi + 1]] + yi + 1]], 8];

                        lerp[lerp[dot2[g00, xa, ya], dot2[g10, xa - 1, ya], u],
                             lerp[dot2[g01, xa, ya - 1], dot2[g11, xa - 1, ya - 1], u], v]],
                        CompilationOptions -> {"InlineExternalDefinitions" -> True}, 
                        CompilationTarget -> "WVM"]];

I had constructed this version of 2D Perlin noise to have a period of $16$ in both of its arguments. Here's how a fundamental piece looks like:

Plot3D[perlin2D[x, y], {x, 0, 16}, {y, 0, 16},
       BoundaryStyle -> None, Mesh -> False, PlotPoints -> 75]

One can use the noise function as is, scale the arguments or the function itself appropriately, or (as is common with how Perlin noise is used) sum so-called "octaves" of them:

Plot3D[perlin2D[x, y] + perlin2D[2 x, 2 y]/2 + perlin2D[4 x, 4 y]/4,
       {x, 0, 16}, {y, 0, 16}, BoundaryStyle -> None, Mesh -> False, PlotPoints -> 75]

Answer 3

Vitaly's answer is correct in that it fantastically produces a splined random disturbace surface. However, I was unable to use it as an initial condition for my NDSolve[...].

Based on whuber's comment and acl's comment, I used :

bsf = Interpolation@Flatten[ Table[{{x, y}, 1 + .05*RandomReal[{-1, 1}]}, {x, 0, L + 1}, {y, 0, L + 1}], 1]

.. as the initial condition for my NDSolve[...]. Yes, NDSolve doesn't like this and complains about a mismatch in the initial and boundary conditions with my favorite ibcinc warning message but I think it is smart enough to reconcile these boundary differences and I am able to solve my partial differential equation satisfactorily.

I tried splining this random surface like whuber and acl suggested but it hasn't worked. If anyone can provide an initial condition that has the SplineClosed->True feature that can be used in NDSolve, that would be truly awesome!

Here is the PDE being solved:

$HistoryLength = 0;
Needs["VectorAnalysis`"]
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Clear[Eq0, EvapThickFilm, h, Bo, \[Epsilon], K1, \[Delta], Bi, m, r]
Eq0[h_, {Bo_, \[Epsilon]_, K1_, \[Delta]_, Bi_, m_, r_}] := \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]h\) + 
        Div[-h^3 Bo Grad[h] + 

      h^3 Grad[Laplacian[h]] + (\[Delta] h^3)/(Bi h + K1)^3 Grad[
        h] + 
            m (h/(K1 + Bi h))^2 Grad[h]] + \[Epsilon]/(
          Bi h + K1) + (r) D[D[(h^2/(K1 + Bi h)), x] h^3, x] == 0;
SetCoordinates[Cartesian[x, y, z]];
EvapThickFilm[Bo_, \[Epsilon]_, K1_, \[Delta]_, Bi_, m_, r_] := 
    Eq0[h[x, y, t], {Bo, \[Epsilon], K1, \[Delta], Bi, m, r}];
TraditionalForm[
    EvapThickFilm[Bo, \[Epsilon], K1, \[Delta], Bi, m, r]];





L = 2*92.389; TMax = 3100*100;
Off[NDSolve::mxsst];
Clear[Kvar];
bsf = Interpolation@
   Flatten[Table[{{x, y}, 1 + .05*RandomReal[{-1, 1}]}, {x, 0, 
      L + 1}, {y, 0, L + 1}], 1];
Kvar[t_] :=  Piecewise[{{1, t <= 1}, {2, t > 1}}]
(*Ktemp = Array[0.001+0.001#^2&,13]*)
hSol = h /. NDSolve[{
          (*Bo,\[Epsilon],K1,\[Delta],Bi,m,r*)

          EvapThickFilm[0.003, 0, 1, 0, 1, 0.025, 0],
          h[0, y, t] == h[L, y, t],
          h[x, 0, t] == h[x, L, t],
          (*h[x,y,0] == 1.1+Cos[x] Sin[2y] *)

          h[x, y, 0] == bsf[x, y]
          },
        h,
        {x, 0, L},
        {y, 0, L},
        {t, 0, TMax},
        Method -> {"BDF", "MaxDifferenceOrder" -> 1},
        MaxStepFraction -> 1/50
        ][[1]]

And the initial condition being plotted: