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.

L = 4*π;
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.


1 Answer 1


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

Therefore $f=∑_{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 π x / lL] + 
                      b[i, j] Sin[2 i π x / lL]) (c[i, j] Cos[2 j π y / lL] + 
                      d[i, j] Sin[2 j π 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 π;
f = fn[5, 5]/5;
Plot3D[Evaluate@f, {x, 0, lL}, {y, 0, lL}]

enter image description here

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

  • $\begingroup$ @ 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. $\endgroup$
    – Enter
    Mar 4, 2015 at 10:07
  • $\begingroup$ 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}] $\endgroup$
    – djp
    Mar 4, 2015 at 13:27
  • $\begingroup$ @ 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?? $\endgroup$
    – Enter
    Mar 4, 2015 at 13:46
  • $\begingroup$ 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. $\endgroup$
    – djp
    Mar 4, 2015 at 22:25
  • $\begingroup$ @ 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. $\endgroup$
    – Enter
    Mar 5, 2015 at 1:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.