Re: FFT in mathematica
- To: mathgroup at smc.vnet.net
- Subject: [mg113688] Re: FFT in mathematica
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Mon, 8 Nov 2010 03:36:58 -0500 (EST)
----- Original Message -----
> From: "Nasser M. Abbasi" <nma at 12000.org>
> To: mathgroup at smc.vnet.net
> Sent: Sunday, November 7, 2010 4:11:29 AM
> Subject: [mg113675] Re: FFT in mathematica
> On 11/6/2010 2:59 AM, Gazi Habiba Akter wrote:
> > Hi,
> > I don't understand how fft algorithm works.
> > but i have to solve a problem by Mathematica code .I am not math and
> > physics student it need my thesis work.
> > it will be ok if any body explain by a simple function.
> > .I tried to solve this problem using fourier transform in
> > mathematica
> > directly but it does not work.I think i need to modify my function
> > or
> > need to know how FFT works for my problem. Here is my code
> > ClearAll["Global'*"]
> > z = 0;
> > n = 100;
> > v1 = 0.5*(1 + Tanh[n*t]);
> > x = (2*t (1 + z))^0.5;
> > v2 = Hypergeometric1F1[.5, 1, -2 x];
> > v3 = Exp[x - t];
> > r1 = v1*v2*v3;
> > f = FourierTransform[r1, t, v];
> > w1 = InverseFourierTransform[f*Exp[z/2*(1 - I*v)], v, t]
> >
> > I need the plot of w1 with respect to t for different value of Z.
> > here
> > r1 is only for z=0.
> > it will be great for me if anybody help me.
> > Thanks,
> > rupa
> >
>
>
> Do you have to do this symbolically? The kernel function you have
> there
> seems to have disconuity at t=0.
>
> Why not simply using numerical solution? using FFT (i.e. Fourier).
>
> Here is one way:
>
> ClearAll["Global'*"]
>
> rupa[z_,t_]:=Module[{n,v1,v2,v3,f,w1,x,r1},
> n=100;
> v1=0.5*(1+Tanh[n*t]);
> x=(2*t (1+z))^0.5;
> v2=Hypergeometric1F1[.5,1,-2 x];
> v3=Exp[x-t];
> r1=Chop[v1*v2*v3];
> f=Fourier[r1];
> w1=InverseFourier[f*Exp[z/2*(1-I*v)]]
> ];
>
> t=Range[-Pi,Pi,Pi/100.]; (*pick discrete values*)
> z=0; (*select z *)
> w1=rupa[z,t];
>
> ListPlot[Inner[List,t,w1,List],
> Frame->True,
> FrameLabel->{{"w1",None},{"t","InverseForurier demo"}}
> ]
>
> --Nasser
This method is much faster than what I'll show below. But it has a few issues. One is that it requires adjustment for correct handling of z!=0. The line
w1=InverseFourier[f*Exp[z/2*(1-I*v)]]
should probably be
w1 = Re[InverseFourier[f*Exp[z/2*(1 - I*t*1/tlength)]]
where tlength is the size of the t range (2*Pi, in Nasser's example). Why take real part? See third issue below.
A second issue involves approximation artifacts (and my approach will also suffer from this, to an extent). Because we use a discrete Fourier transform, we are (implicitly) periodicizing the function. In consequence, we cannot expect sensible results when we check near the boundaries (+-Pi, in the case above).
A third issue, related to the second, is we now can get imaginary parts in the artifacts. When I make the range {-4*Pi,4*Pi} these seem to be reduced e.g. at z=1. I'm guessing this is some flavor of Gibb's phenomenon, behind the scenes.
The code below does different numerica approximating, using quadrature. I also discard v1[t] in some places (where I try to set up an exact approximation to a subproblem). The reason here is that it is quite close to a simple UnitStep, so I approximate by 1 for t>=0 and 0 otherwise (as reflected in bounds of integration).
A check of the series expansions will show that this function drops rapidly to zero, and restricting integration range to 10 seems a useful cutoff.
n = 100;
v1[t_] := (1 + Tanh[n*t])/2;
x[t_, z_] := (2*t (1 + z))^(1/2);
v2[t_, z_] := Hypergeometric1F1[1/2, 1, -2 x[t, z]];
v3[t_, z_] := Exp[x[t, z] - t];
r0[t_, z_] := v2[t, z]*v3[t, z];
r1[t_, z_] := v1[t]*r0[t, z]
r2[t_, z_] = Normal[Series[r0[t, z], {t, 0, 30}]];
I tried three approximations, as below. The second one is fastest to evaluate, but has horrible numeric conditioning (needs very high precision). I think this may be an issue of cancellation error in numerical evaluation of a series expansion, but did not really check too carefully. There may well be a smarter way that avoids the trouble and retains the speed.
ff[v_?NumericQ, z_?NumericQ] :=
NIntegrate[r1[t, z]*Exp[I*v*t], {t, -20, 20}]/Sqrt[2*Pi]
ff2[v_, z_] = Integrate[r2[t, z]*Exp[I*v*t], {t, 0, 10}]/Sqrt[2*Pi];
ff3[v_?NumericQ, z_?NumericQ] :=
NIntegrate[r2[t, z]*Exp[I*v*t], {t, 0, 10}]/Sqrt[2*Pi];
w1[t_?NumericQ, z_?NumericQ] :=
NIntegrate[ff[v, z]*Exp[z/2*(1 - I*v)]*Exp[-I*t*v], {v, -10, 10}]/
Sqrt[2*Pi]
w2[t_?NumericQ, z_?NumericQ] :=
NIntegrate[ff2[v, z]*Exp[z/2*(1 - I*v)]*Exp[-I*t*v], {v, -10, 10},
WorkingPrecision -> 200, PrecisionGoal -> 40, AccuracyGoal -> 40]/
Sqrt[2*Pi]
w3[t_?NumericQ, z_?NumericQ] :=
NIntegrate[ff3[v, z]*Exp[z/2*(1 - I*v)]*Exp[-I*t*v], {v, -10, 10}]/
Sqrt[2*Pi]
Some evaluations:
In[194]:= w3[1/2, 1]
Out[194]= 1.42882 - 1.72726*10^-13 I
In[217]:= w3[3, 1]
Out[217]= 0.464579 + 0. I
w1 seems to agree in the evaluation below, quadrature messages notwithstanding.
In[240]:= Timing[w1[1/2, 1]]
During evaluation of In[240]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>
During evaluation of In[240]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in v near {v} = {0.429377}. NIntegrate obtained 3.5793+0. I and 0.00001583496868014035` for the integral and error estimates. >>
Out[240]= {35.677, 1.42793 + 0. I}
I do not recommend using my code as the primary means for what you are doing. But it might be useful as a means of testing modifications of Nasser Abbasi's (hugely faster) code, in order to get some confidence as to what is correct, what ranges are "safe" from the periodicization effect, etc. Possibly Abbasi or another respondent will be able to give a more definitive account of what is the correct way to use the discrete transform to approximate the continuous one. I'm just guessing, above. Might be right, might be wrong, might be Green Eggs and Ham.
Daniel Lichtblau
Wolfram Research