# Small random disturbance of a flat surface

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$

$HistoryLength = 0; Needs["VectorAnalysis"] Needs["DifferentialEquationsInterpolatingFunctionAnatomy"]; 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 ][]  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. >> • What is the matter with the "white noise" Plot3D[1 + RandomReal[{-0.05, 0.05}], {x, 0, 1}, {y, 0, 1}]? That would be valid as an initial condition provided it were repeatable, which can be accomplished by memoizing it: f[x_, y_] := f[x, y] = RandomReal[]; Plot3D[f[x, y], {x, 0, 1}, {y, 0, 1}]. Just watch out for uncontrolled growth in RAM used by f! Dec 4, 2012 at 22:37 • @whuber Good point. However, this seems rather computationally intensive.. like you point out. Dec 4, 2012 at 22:41 • Yes, but if that's what's intended... . I rather suspect, though, that you might want to refine your concept of a "random disturbance." You are asking for a random spatial field and those have structure. What structure are you looking for? I have described one that has a particularly simple covariance function, but (consequently) the realizations are not even continuous. Dec 4, 2012 at 22:45 • @whuber Sorry, but I barely understand the tech language that you just used. I've included an example with the bspline thingy and it errors out. Dec 4, 2012 at 22:46 • @drN I don't know but if you do bsf = Interpolation@Flatten[ Table[{{x, y}, 1 + .05*RandomReal[{-1, 1}]}, {x, 0, L + 1}, {y, 0, L + 1}], 1]; and then have as an ic h[x, y, 0] == bsf[x, y] it works (except I didn't fix the boundary conditions correctly) – acl Dec 4, 2012 at 23:49 ## 3 Answers 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}] • So is this correct: RandomReal creates a random number array in the range (0,1) which is 30x30. Then this is fed into BSplineFunction to create a spline curve? Dec 4, 2012 at 22:34 • @drN Yes, but spline surface (not curve) which you cna use as analitic function: find derivative or feed into a NDSolve. Dec 4, 2012 at 22:37 • Well, using this as initial condition didn't give me agreeable results. My initial condition is h[x,y,0]==BSplineFunction[RandomReal[1, {10, 10, 1}]] for quite a complicated PDE... I'll see if I can include a simple working problem. Dec 4, 2012 at 22:40 • @drN Do you have to take care of boundary conditions? Dec 4, 2012 at 22:41 • @ VitaliyKaurov, I try to use your answer f = BSplineFunction[RandomReal[1, {30, 30, 1}], SplineClosed -> True] to produce a small random disturbance of a plane in a square domain of [0,4*pi]*[0,4*pi] which has periodic boundary condition. But when I try, say f[9,0]==f[9,4*pi], MMA gives me False. I know you produce a disturbed surface in the domain [0,1]*[0,1], but how to get such surface in the domain, say [0,4*pi]*[0,4*pi]? Mar 2, 2015 at 12:20 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] 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["DifferentialEquationsInterpolatingFunctionAnatomy"];
Clear[Eq0, EvapThickFilm, h, Bo, \[Epsilon], K1, \[Delta], Bi, m, r]
Eq0[h_, {Bo_, \[Epsilon]_, K1_, \[Delta]_, Bi_, m_, r_}] := \!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]h\) +

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}];
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
][]
`

## And the initial condition being plotted: 