# Generate a small random disturbance on a flat surface for NDSlove

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 between 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.

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 and Vitaliy's answer, 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 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 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.

Every function of the following form satisfies your boundary conditions: $C_{i,j}=\left(a_{i,j}\cos(2i\pi x)+b_{i,j}sin(2i\pi x\right)\left(c_{i,j}\cos(2j\pi y)+d_{i,j}sin(2j\pi y\right)$

Therefore $f=\sum_{i,j} C_{i,j}$ also satisfies, for any random coefficients $a$, $b$, $c$, $d$. Here's a possible implementation:

component[i_, j_] := (a[i, j] Cos[2 i \[Pi] x / lL] +
b[i, j] Sin[2 i \[Pi] x / lL]) (c[i, j] Cos[2 j \[Pi] y / lL] +
d[i, j] Sin[2 j \[Pi] y / lL])
fn[imax_, jmax_] := Sum[component[i, j], {i, 0, imax}, {j, 0, jmax}]
a[i_, j_] := RandomReal[{-1, 1}]
b[i_, j_] := RandomReal[{-1, 1}]
c[i_, j_] := RandomReal[{-1, 1}]
d[i_, j_] := RandomReal[{-1, 1}]
lL=4 \[Pi];
f = fn[5, 5]/5;
Plot3D[Evaluate@f, {x, 0, lL}, {y, 0, lL}]


It works with your code, as long as you set the initial condition to be h[x,y,0]==f.

• @ djp, thanks for your answer. I have no a fully understand as to f = fn[5, 5]/5. I guess f = fn[5, 5]/6 might be more suitable. Let me try ur answer with my full problem. – Enter Mar 4 '15 at 10:07
• Keep in mind my random coefficients were entirely ad hoc. If you use them I suggest you do something like a[i_,j_] := Exp[-(i+j) RandomReal[{-1,1}] – djp Mar 4 '15 at 13:27
• @ djp, thank for you kindly help! What is your mean by ad hoc? So is this correct: by a[i_,j_] := Exp[-(i+j) RandomReal[{-1,1}], we can make the random coefficients more analytical, which can be fed into NDSlove more suitable or more efficiently?? – Enter Mar 4 '15 at 13:46
• Well I don't know your problem, but it seems likely to me that you want low frequency perturbations (i,j ~ 0, 1, 2) to be more intense than high frequency perturbations. – djp Mar 4 '15 at 22:25
• @ djp, There seems to be a right-square bracket before RandomReal in your earlier comment? The simulation is running. I will give you a feedback in two days. – Enter Mar 5 '15 at 1:35