Re: envelope of an oscillatory InterpolatingFunction
- To: mathgroup at smc.vnet.net
- Subject: [mg49812] Re: envelope of an oscillatory InterpolatingFunction
- From: Paul Abbott <paul at physics.uwa.edu.au>
- Date: Sun, 1 Aug 2004 18:48:44 -0400 (EDT)
- Organization: The University of Western Australia
- References: <cefh67$bs5$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
In article <cefh67$bs5$1 at smc.vnet.net>,
"Stergios J. Papadakis" <stergios.papadakis at jhuapl.edu> wrote:
> I use NDSolve and get an oscillatory
> InterpolatingFunction with a slowly
> varying envolope. I would like to
> extract just the envelope. Right
> now I am FFTing the InterpolatingFunction,
> using NMaximize to find the peak and
> therefore the carrier frequency.
> Then I step through the InterpolatingFunction
> using a While loop, one period at a time,
> and use NMaximize to find the maximum on
> that each period. This is really slow, and sometimes
> it just obviously gets it wrong and I can't
> really figure out why.
>
> Is there a better way that I just don't know about?
Why not compute the values of the extrema of the InterpolatingFunction
and then use FindFit to determine the envelope? For example,
nsol = First[NDSolve[{10 y''[x]+y'[x]+5 y[x]==0,y[0]==1,y'[0]==0},
y, {x, 0, 50}]]
Plot[Evaluate[y[x] /. nsol], {x, 0, 50}, PlotRange -> All]
Here is some code for finding the roots of a function over a specified
range:
Needs["Utilities`FilterOptions`"]
RootsInRange[fn_, {x_, xmin_, xmax_}, opts___] :=
Module[{p, pts, x, f = Function[x, Evaluate[fn]]},
p = Plot[f[x], {x, xmin, xmax}, Compiled -> False,
Evaluate[FilterOptions[Plot, opts]]];
pts = Cases[First[p], Line[{x__}] -> x, Infinity];
pts = Map[First, Select[Split[pts,
Sign[Last[#2]]==-Sign[Last[#1]]&],Length[#1]==2& ],{2}];
(FindRoot[f[x] == 0, {x, Sequence @@ ##1},
Evaluate[FilterOptions[ FindRoot, opts]]]&) /@ pts]
Compute the (absolute) values of the extrema of the
InterpolatingFunction:
{x, Abs[y[x]]} /. RootsInRange[Evaluate[y'[x] /. nsol],{x,0,50}] /. nsol
and then fit to find the Envelope (here assuming exponential form):
FindFit[%, a Exp[-b x] + c, {a, b, c}, x]
Cheers,
Paul
--
Paul Abbott Phone: +61 8 9380 2734
School of Physics, M013 Fax: +61 8 9380 1014
The University of Western Australia (CRICOS Provider No 00126G)
35 Stirling Highway
Crawley WA 6009 mailto:paul at physics.uwa.edu.au
AUSTRALIA http://physics.uwa.edu.au/~paul