Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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 \[Pi] x/L] + Sin[2 \[Pi] x/L]) Cos[
    2 \[Pi] 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.

Working example:

$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];
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] == 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 BSpline 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. >>

share|improve this question
    
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! –  whuber Dec 4 '12 at 22:37
    
@whuber Good point. However, this seems rather computationally intensive.. like you point out. –  drN Dec 4 '12 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. –  whuber Dec 4 '12 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. –  drN Dec 4 '12 at 22:46
1  
@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 '12 at 23:49

3 Answers 3

up vote 10 down vote accepted

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]]

enter image description here

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}]

enter image description here

share|improve this answer
    
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? –  drN Dec 4 '12 at 22:34
    
@drN Yes, but spline surface (not curve) which you cna use as analitic function: find derivative or feed into a NDSolve. –  Vitaliy Kaurov Dec 4 '12 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. –  drN Dec 4 '12 at 22:40
    
@drN Do you have to take care of boundary conditions? –  Vitaliy Kaurov Dec 4 '12 at 22:41
    
It probably is as I need to use Periodic boundary conditions. However, mathematica didn't scream at me that the initial and bound conditions are inconsistent like it would usually do for bad ics/ –  drN Dec 4 '12 at 22:42

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]

basic Perlin noise

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]

three octaves of Perlin noise

share|improve this answer

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:

enter image description here

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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