Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Finding a sine wave

  • To: mathgroup at smc.vnet.net
  • Subject: [mg95390] Re: [mg95341] Finding a sine wave
  • From: DrMajorBob <btreat1 at austin.rr.com>
  • Date: Sat, 17 Jan 2009 05:30:11 -0500 (EST)
  • References: <200901161107.GAA14041@smc.vnet.net>
  • Reply-to: drmajorbob at longhorns.com

Your constraint on the wavelength will NOT make solutions unique, I think,  
and not all data will yield a solution at all.

Here are three solvers you might try; the first is yours, except that I've  
put constraints on the wavelength. Your version set a starting point,  
only, for the wavelength.

Clear[findSin, fitSin, minSin]
findSin[xx_List, yy_List] :=
  findSin[xx, yy] =
   Block[{a, ph, wl, x},
    a Sin[ph + 2 Pi x/wl] /.
     FindRoot[
      yy - a Sin[ph + 2 Pi xx/wl] //
       Evaluate, {{a, Max[Abs[yy]]}, {wl, 2 Max@xx, Max@xx,
        10 Max@xx}, {ph, 0}}]]
fitSin[xx_List, yy_List] := fitSin[xx, yy] = Block[{a, ph, wl, x},
    a Sin[ph + 2 Pi x/wl] /.
     FindFit[Rest@Transpose@{xx, yy}, {a Sin[ph + 2 Pi x/wl],
       wl > Max@Abs@xx, 0 < 2 Pi ph < Max@Abs@xx}, {a, ph, wl}, x]
    ]
minSin[xx_List, yy_List] := minSin[xx, yy] = Block[{a, ph, wl, x},
    a Sin[ph + 2 Pi x/wl] /.
     Last at NMinimize[{#.# &[Thread[yy - a Sin[ph + 2 Pi xx/wl]]],
        wl > Max@xx, 0 < ph < (1/2 Pi) Max@xx}, {a, ph, wl}]
    ]

xx=Sort[Join[{0},RandomReal[{0,10},2]]];
yy=RandomReal[{-25,25},3];
f[x_]=findSin[xx,yy]
g[x_]=fitSin[xx,yy]
h[x_]=minSin[xx,yy]

-27.4879 Sin[0.456353-1.55521 x]
FindFit::eit: The algorithm does not converge to the tolerance of  
4.806217383937354`*^-6 in 500 iterations. The best estimated solution,  
with feasibility residual, KKT residual or complementary residual of  
{1.0772*10^-10,0.0000226396,3.80526*10^-11}, is returned. >>
-37.2328 Sin[0.527053+0.0356086 x]
-27.4879 Sin[2.68524+1.55521 x]

Plot[{f@x, g@x, h@x}, {x, 0, 10},
  Epilog -> {PointSize[0.02], Point /@ Transpose@{xx, yy}},
  PlotRange -> {{0, 10}, {-25, 25}}]

minSin will (almost always?) return a solution without errors, but it  
won't always fit the data points. The other two will fail in ALL SORTS of  
ways.

In only three random draws, I found an example for which none of the  
solvers actually fit the data:

xx={0, 1.5565, 6.32643}
yy={0.305144, 3.21477, 17.9622}
f[x_] = findSin[xx, yy]
g[x_] = fitSin[xx, yy]
h[x_] = minSin[xx, yy]

FindRoot::reged: The point {21.1716,63.2643,0.000205665} is at the edge of  
the search region {6.32643,63.2643} in coordinate 2 and the computed  
search direction points outside the region. >>
21.1716 Sin[0.000205665+0.0993164 x]
16.5886 Sin[1.00688+0.993164 x]
-9.61001 Sin[4.39622+0.993164 x]

findSin fit the first two points, fitSin fit NONE of the points, and  
minSin fit only the second point.

I'm pretty sure NO Sin function fits those data points, with a wavelength  
longer than the interval.

Bobby

On Fri, 16 Jan 2009 05:07:49 -0600, Hugh Goyder  
<h.g.d.goyder at cranfield.ac.uk> wrote:

> An experiment will give me three coordinates which lie on a sine wave.
> I have to find the sine wave efficiently. There are three unknowns the
> sine wave amplitude, A, the wavelength of the sine wave, L,  and the
> phase, ph. I also know that the wavelength is larger than the interval
> containing my measurement points. I think this condition removes
> possible multiple solutions.
>
>  Below I give two methods which need improving. The first method uses
> FindRoot. This methods works about 60% of the time. I give six
> examples where it fails. The failure may be due to poor initial
> guesses. In the second method I use FindFit. This is not quite the
> correct method because I have an equal number of equations and
> unknowns. Some of the failures here, I think, are due to there being
> no error for the algorithm to work with. I give examples of failures.
> I have also tried FindInstance, Reduce and NSolve but I don't think
> these are  appropriate.
>
> Here are some questions
>
> 1. I would really like a symbolic solution rather than an iterative
> one. Is such a solution possible?
> 2. Can anyone improve on the methods below to make them more robust?
> 3. I have some control over my x-locations. How can I work out best x
> locations given an estimate of the wavelength L?
>
> Many thanks for all answers.
>
>
> xx = Sort[Join[{0}, RandomReal[{0, 10}, 2]]]
>
> yy = RandomReal[{-25, 25}, 3]
>
> sol = FindRoot[{yy[[1]] - A*Sin[ph],
>    yy[[2]] - A*Sin[2*Pi*(xx[[2]]/L) + ph],
>        yy[[3]] - A*Sin[2*Pi*(xx[[3]]/L) + ph]}, {{A,
>     Max[Abs[yy]]}, {L, Max[xx]}, {ph, 0}}]
>
> Plot[Evaluate[A*Sin[2*Pi*(x/L) + ph] /. sol], {x, 0, 10},
>    Epilog -> {PointSize[0.02], (Point[#1] & ) /@
>     Transpose[{xx, yy}]}, PlotRange -> {{0, 10}, {-25, 25}}]
>
> FindRootFailures = {{{0, 3.781264462608982,
>      3.797055100562352}, {-22.948348087068737, 1.4581744078038472,
>            -14.242676740704574}}, {{0, 4.424069570670131,
>      8.861716396743098}, {-6.444538773775843, -10.787309362608688,
>            15.579334094330942}}, {{0, 4.424069570670131,
>      8.861716396743098}, {-5.944536471563289, 10.627536107497393,
>            -21.294232497202316}}, {{0, 1.3680556047878967,
>      8.546115267250002}, {4.128528863849845, 2.9017293848933923,
>            -22.51610539815371}}, {{0, 1.3738869718371616,
>      5.689309423462079}, {20.553993773523437, -16.972841620064592,
>            16.61185061676568}}, {{0, 1.0408828831133632,
>      8.484645515699821}, {-3.7267589861478045, 4.1016610850387,
>            -24.600635061804443}}};
>
> (({xx, yy} = #1;
>     sol = FindRoot[{yy[[1]] - A*Sin[ph],
>        yy[[2]] - A*Sin[2*Pi*(xx[[2]]/L) + ph],
>
>        yy[[3]] - A*Sin[2*Pi*(xx[[3]]/L) + ph]}, {{A,
>         Max[Abs[yy]]}, {L, Max[xx]}, {ph, 0}}];
>         Plot[Evaluate[A*Sin[2*Pi*(x/L) + ph] /. sol], {x, 0, 10},
>
>      Epilog -> {PointSize[0.02], (Point[#1] & ) /@
>         Transpose[{xx, yy}]},
>      PlotRange -> {{0, 10}, {-25, 25}}]) & ) /@
>    FindRootFailures
>
> xx = Sort[Join[{0}, RandomReal[{0, 10}, 2]]]
>
> yy = RandomReal[{-25, 25}, 3]
>
> sol = FindFit[
>   Transpose[{xx, yy}], {A*Sin[2*Pi*(x/L) + ph],
>    Max[Abs[xx]] < L}, {A, L, ph}, x]
>
> Plot[Evaluate[A*Sin[2*Pi*(x/L) + ph] /. sol], {x, 0, 10},
>    Epilog -> {PointSize[0.02], (Point[#1] & ) /@
>     Transpose[{xx, yy}]}, PlotRange -> {{0, 10}, {-25, 25}}]
>
> FindFitFailures = {{{0, 2.463263668380835,
>      4.3190892163093615}, {-1.8407827676541144, -8.736574079785198,
>            12.661520984622932}}, {{0, 3.894521446091823,
>      9.937403619870642}, {12.381369822165155, 15.840399165432128,
>            -6.914634137727327}}, {{0, 8.087369725271945,
>      8.343899312282815}, {24.73795103895976, -13.396248713970305,
>            7.150470311065216}}, {{0, 1.5480031866178834,
>      9.575255260205617}, {-8.163720246784278, 13.373468882892958,
>            -18.018462091502098}}, {{0, 0.6152784485896601,
>      0.818296772602134}, {-1.858698140836046, 3.695113783904491,
>            -1.3186989232026658}}, {{0, 4.743093154057316,
>      6.028314583406327}, {-16.60893277597446, -17.413392198343093,
>            -18.54081837965986}}};
>
> (({xx, yy} = #1;
>     sol = FindFit[
>       Transpose[{xx, yy}], {A*Sin[2*Pi*(x/L) + ph],
>        Max[Abs[xx]] < L}, {A, L, ph}, x];
>         Plot[Evaluate[A*Sin[2*Pi*(x/L) + ph] /. sol], {x, 0, 10},
>
>      Epilog -> {PointSize[0.02], (Point[#1] & ) /@
>         Transpose[{xx, yy}]},
>      PlotRange -> {{0, 10}, {-25, 25}}]) & ) /@
>    FindFitFailures
>
>
>
>



-- 
DrMajorBob at longhorns.com


  • Prev by Date: Re: Solve / NSolve
  • Next by Date: Re: Mathematica and Access
  • Previous by thread: Finding a sine wave
  • Next by thread: Re: Finding a sine wave