Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

This question already has an answer here:

I am solving a differential equation numerically and the output is an oscillating function with the amplitude of the oscillation decaying in time. I would like to extract the power law governing this amplitude.

Let's make this concrete. Let's say I am solving the system

x''[t] + 2t^(-1) x'[t] + t^(-2) x[t] == 0

which has solutions

x= C[1] t^(-1/2) Cos[Sqrt[3]/2 Log[t]]

and the corresponding $\sin$. I would like to extract the power law $t^{-1/2}$.

Of course, my system is much more complicated and I have to solve it numerically. I get as my solution an interpolating function. I was thinking of generating a table of function values, sorting according to the maxima and then fitting a power law, but maybe there is a better way?


share|improve this question

marked as duplicate by Jens, m_goldberg, Kuba, Sjoerd C. de Vries, Artes Aug 6 '13 at 10:26

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

I think your own suggestion is a good one! – Ali Jun 27 '13 at 10:09
You may find this discussion relevant. – Leonid Shifrin Jun 27 '13 at 10:27
If you're dealing entirely with cases that have analytic solutions, you can likely find the envelope analytically as well: (Just noticed the last part of the post.. oops. Still worth mentioning, though) – Corey Kelly Jun 27 '13 at 11:19
There is also an enveloping function to be found in this question:… – Sjoerd C. de Vries Aug 6 '13 at 6:17

Here is a general envelope method developed for this SystemModeler industry example:

FunctionEnvelope[f_, {t_, a_, b_}, n_: 40] := 
 Module[{seeds, x, y, points, progress = 0, tempf, union}, 
  seeds = Rescale[Range[0, 1, 1/n] + 1/(2 n), {0, 1}, {a, b}]; 
  points = Last[
        Quiet@Check[progress++; {y, x} = FindMaximum[Abs[f], {t, x0}];
           x = t /. x; If[a <= x <= b, Sow[{x, y}]], Null], {x0, 
         seeds}], ProgressIndicator[progress, {0, Length[seeds]}]]]]];
  union[] :=
   points = 
     SameTest -> (Abs[First[#1] - First[#2]]/
         Replace[Max[Abs[First[#1]], Abs[First[#2]]], 
          u_ /; u == 0 :> 1] < 10^-6 &)];
  tempf = Interpolation[points];
  points = Quiet[Join[{{a, tempf[a]}}, points, {{b, tempf[b]}}]];

It's used like this:

f[x_] := x^2 Sin[x] Sin[x^2]

g = FunctionEnvelope[f[x], {x, 0, 15}, 100];

Now g[x] is the automatically determined envelope function:

Plot[{f[x], g[x], -g[x]}, {x, 0, 15}]

f[x] and automatically determined envelope g[x]

This works by smoothly interpolating between the local maxima of the function.

share|improve this answer

The idea of an envelope is often quite clear, but it is a slippery thing to try and define it mathematically. In audio, the envelope is often used to help locate start and ending points for events, and is often associated with two parameters: a rise time (how fast the envelope is allowed to grow) and a decay time (how fast the envelope is allowed to die away. Concretely, one can calculate the envelope as a combination of two linear systems with two different time constants.

decay = 0.006; rise = 0.2;
filt[z_, u_] := Max[decay z + (1 - decay) u, rise z + (1 - rise) u];
env = FoldList[filt, 0, Abs[data]];

The filt function implements the two linear systems and combines them with Max. The data is a list (in audio it might be a single channel extracted from a wav file), and the two parameters (both must be between zero and one) specify the rise and decay times.

Here is this enveloping operation in action. The light red is the data and the dark blue is the calculated envelope.

enter image description here

share|improve this answer

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