Message Boards Message Boards

Index an arbitrary array in NDSolve to simulate white noise?

Posted 6 years ago

I'm trying to introduce a randomized element in a simulation, using NDSolve to do the integration. The random element is supposed to simulate a random thermal force. The way I've tried to implement this is by creating an array of random values and then trying to select values from the array during each time step. The following is a snippet of my current code.

(*Initialising some time values that will allow me to select particular values from my random array*)
Tfinal = 10;
SampleFreq = 2500;
NoT = Tfinal*SampleFreq + 1;
dt = Tfinal/(NoT - 1);

(*Array of random values and the function I'll be using in NDSolve*)
WN = RandomVariate[NormalDistribution[0, 1], NoT];
WNoise[t_] := Indexed[WN, Floor[Floor[t/dt]] + 1]

(*NDSolve function using white noise function*)
sol = x[t] /. 
  NDSolve[{D[x[t], t, t] == -(w^2) x[t] + WNoise[t], x[0] == 0, 
    x'[0] == 2 }, x, {t, 0, 10}]

The function above (without white noise) is the classic F = -kx equation. I've just tried to introduce a white noise driving force force to it as well. Unfortunately, it doesn't work when the WNoise function is added. The error message I typically get is the following.

Indexed::ind: The index 1/2 is not a nonzero integer.
NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`.

I'm just wondering if anyone has faced a similar problem and if this method of introducing white noise is possible or if I'm just being silly with this implementation.

Thank you for your help.

POSTED BY: Sandip Roy
7 Replies
Anonymous User
Anonymous User
Posted 6 years ago

I forgot the most important principle (always a bad thing to do).

Observe F sin(wt+B) as the (counter/added) force function to -kx in time, and that it must have some ENERGY (work) (or F is 0, nothing).

If the force is limited, the effect of it will be limited if there is any DAMPING (imagine a spring oscillating in water, the water will stop the motion of the oscillating spring. If there were a "added force" to the oscillation it would need to BE MORE THAN the water damping to have any OVERALL IMPACT).

The overall analysis should be done first, which is that the energy of the random "addition" matters, and more so if the damping is greater than randomness matters more than that.

POSTED BY: Anonymous User
Anonymous User
Anonymous User
Posted 6 years ago

introduce a randomized element in a simulation base simulation is of F = -kx, with heat introducing a random element

it's vague to me what you are simulating. i see frequency and unsure what it is the frequency of (thus, how heat could physically effect it). is it a spring? is random heat lengthening/softening the spring but random cooling also occurs? (but that situation doesn't make sense to simulate)

but if the randomness were "average random", then overall there would be no change in derivatives or integrals

without knowing what you are trying to simulate suggesting answers doesn't make sense to me

since both eq. above are solved in many books, there would be no need to solve them (you could look them up). but mm can do both forms I'm sure Neil Singer (1st responder) knows better than me as to answering.

I see your equation as form (t->x, using y,x as typical, and Q(x) representing the random (equation of in x if we call it that, and what kind matters)):

 y' + c y = Q(x)

that's certainly solvable if entered into DSolve correctly. and you can solve for Yp for each Q(x) or for fun entertain a family of Q(x) as Yp.

however for y = -kx, the "standard" eq. (time, motion, position) to use is:

y'' + w0^2 y = F sin(wt + B) (* forced undamped motion is the example, w is a constant *)

which is stolen from p. 365 of Ord. Diff. Eq. (tenenbaum, pollard). you might then wonder at the solutions of that if B is random. but it would be best to not do that and speculate for that particular case that it would result in random phase change

(if you didn't know, Y(x)=Yc+Yp, where Yp =0 if Q(x)=0, that is solutions for Q(x) are simply added to the general solution Yc which is when Q(x)=0)

MAYBE this will help answer the question. for y''+w0^2y=F sin(wt+B), if w0 is close to w, wild surges occur, otherwise we may see damping (a ring, peters out), or sin some situations a gradual rise or fall in a range (not a full complete sine wave of motion).

SUMMARY of MAYBE: if you take the opposite of what Neil Singer said (force the randomness to have a frequency) and do NOT force a frequency, then you can say nothing of the result. it could be zero by chance, increase or decrease to zero, it could decrease to zero then revive itself.

POSTED BY: Anonymous User

Try searching http://reference.wolfram.com/language/ref/ItoProcess.html (the doc. page for ItoProcess) for "white" (or "white noise"). You'll should find an example very similar to what you're trying to do.

POSTED BY: Michael Rogers

Sandip,

Glad I could help.

You created a list of random samples and you wanted them at SampleFreq spacing to represent noise with a certain frequency content. To do this I had to create a continuous time function with your noise. Interpolation[] takes a list of pairs in the form of {time, value}. I used a "trick" to create the list of pairs:

  1. I created the time base as a vector (list) of times with Range[0,10,1/SampleFreq] -- {0, 1/2500, 1/1250,...}
  2. I made a matrix out of it by putting the times and values in a list: {{0, 1/2500, 1/1250,...}, {.111, .222, .333, ...}} This is a matrix with two rows and many columns, However the format is wrong for Interpolate -- I need pairs of points which consist of many rows of two Columns.
  3. I use Transpose[] on the matrix to make it have the correct format. as an example:

    In[1]:= lst = {{1, 2, 3, 4, 5}, {a, b, c, d, e}}
    
    Out[1]= {{1, 2, 3, 4, 5}, {a, b, c, d, e}}
    
    In[2]:= Transpose[lst]
    
    Out[2]= {{1, a}, {2, b}, {3, c}, {4, d}, {5, e}}
    

In matrixform:

enter image description here

I used Interpolation to create a function of time that the integrator can use as any other continuous function in MMA.

Now you may want to consider some options to Interpolation because the default is 3rd order interpolation. If you are modeling a digital system with noise you may want to use InterpolationOrder->0 to get a sample-and-hold version.

The Transpose[] trick is worth remembering because it is frequent that you have lists in the wrong format and need to rearrange things. I call it a "trick" because most people do not think of matrix operations when handling lists but they can be useful.

Regards,

Neil

POSTED BY: Neil Singer

Sandip,

You can't mix a digitized noise signal with a continuous integration that way. The derivatives are evaluated along continuous time points and are not on fixed intervals. You have several options. You can use WhenEvent to trigger an event at each sample point and add in the noise or you can use Interpolation to make your noise continuous. I'll show the later. You also need to decide the frequency of the noise (SampleFreq) because a very high rate will significantly slow down the integration.

Tfinal = 10;
SampleFreq = 2500;
NoT = Tfinal*SampleFreq + 1;
dt = Tfinal/(NoT - 1);

(*NDSolve function using white noise function*)

WN = RandomVariate[NormalDistribution[0, 1], NoT];
noise = Interpolation[Transpose[{Range[0, 10, 1/SampleFreq], WN}]];
sol = x[t] /. 
   NDSolve[{D[x[t], t, t] == -(w^2) x[t] + 100*noise[t], x[0] == 0, 
     x'[0] == 2}, x, {t, 0, 10}];
Plot[sol[[1]], {t, 0, 10}]

Regards,

Neil

POSTED BY: Neil Singer
Posted 6 years ago

Hi, Neil.

Thank you so much! Your solution is exactly what I was looking for.

However, I am quite new to Mathematica and I was wondering if you could explain to me what exactly the function Interpolation[Transpose[{Range[0, 10, 1/SampleFreq], WN}]] is doing. Naively, to me, it looks like it's taking the WN array and interpolating between the different values, with a frequency of SampleFreq from t=0 to t=10. However, I don't quite understand the purpose of the Transpose function.

Thank you again for your help, Neil.

Best Regards,

Sandip

POSTED BY: Sandip Roy

Other responses are probably a better way to go but I want to mention that this can be made to work with very simple modifications.

(1) Give w a value.

(2) Define WNoise so that it only evaluates for explicit numeric values of t.

WNoise[t_?NumberQ] := Indexed[WN, Floor[t/dt] + 1]

(3) Also it would be good to maintain consistency with time parameters, so run NDSolve from 0 to Tfinal.

A possible advantage to this approach is that it does not try to create a differentiable function from the random impulses. Which is no guarantee NDSolve won't treat it as differentiable using numeric finite differences (though I do not think it will do that).

POSTED BY: Daniel Lichtblau
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract