Re: Why this function does not return a single value
- To: mathgroup at smc.vnet.net
- Subject: [mg60180] Re: Why this function does not return a single value
- From: Paul Abbott <paul at physics.uwa.edu.au>
- Date: Tue, 6 Sep 2005 01:26:50 -0400 (EDT)
- Organization: The University of Western Australia
- References: <dfiv9c$p99$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
In article <dfiv9c$p99$1 at smc.vnet.net>,
Marek Bromberek <marek at chopin.physics.mun.ca> wrote:
> I was wondering if somebody could tell me what am I doing wrong here.
A number of things:
[1] The syntax is incorrect. You cannot use (a[1]*_)?NumericQ. This
means a[1] times a blank, which is not what you intend.
If you enter
(x_)?NumericQ // FullForm
and
(a[1]*_)?NumericQ // FullForm
you will see the difference.
[2] A much better syntax is to define a test for numerical lists, say
num[c_List] := And @@ (NumericQ /@ c)
and then use it in your definition, say
VoigtSum[x_?NumericQ, a_List?num, b_List?num, ...] := ...
> I am constructing a function (for fitting purposes) which is a sum of 15
> Voigt functions + one Gaussian peak and and exponential background. However
> when I want to check if that function works and try to evaluate it it does
> not return a value.
[3] There is no need to restrict attention to 15 Voight functions or to
test the length of each list. Now your arguments are lists, instead
write your definition so that it works for _any_ number of Voight
functions (mapping the function over the lists of arguments), and use
the last two elements of the lists a, b, and \[Delta]G to give the
Gaussian peak and exponential background, i.e.,
a[[-2]] Exp[(-Log[2]) ((x - b[[-2]])^2/\[Delta]G[[-2]]^2)] +
a[[-1]] Exp[-x/b[[-1]]]
[4] You have
> Sum[a[i]*((2.*Log[2.]*\[Delta]L[i])/(N[Pi^(3/2)]*\[Delta]G[i]))*
> NIntegrate[Exp[-t^2]/((Sqrt[Log[2.]]*(\[Delta]L[i]/\[Delta]G[i]))^2 +
> (Sqrt[4.*Log[2.]]*((x - b[i])/\[Delta]G[i]) - t)^2), {t, -Infinity,
> Infinity}], {i, 1, 15}] +
> a[16]*Exp[(-Log[2.])*((x - b[16])^2/\[Delta]G[16]^2)] +
> a[17]*Exp[-x/b[17]]
There is no need -- and it is usually a bad idea -- to numericalize
constants. Since you are computing a numerical integral, all exact
numerical values will be coerced to numerical ones (of the same
precision). Try entering
2. Pi
to see what I mean.
[5] Although Mathmatica cannot compute your integral directly,
Integrate[Exp[-t^2]/((Sqrt[Log[2.]] (\[Delta]L[i]/\[Delta]G[i]))^2 +
(Sqrt[4.*Log[2.]]*((x - b[i])/\[Delta]G[i]) - t)^2),
{t, -Infinity, Infinity}]
or the equivalent integral
Integrate[Exp[-t^2]/(a^2 + (b - t)^2), {t, -Infinity, Infinity}]
it is possible to compute this in closed form (at least for a and b
real). I get
int[a_, b_] = -((1/(2 a)) (I (((-I) Pi Erf[a + I b] -
Log[1/(I a - b)] + Log[1/(b - I a)])/E^(b - I a)^2 -
(I Pi Erf[a - I b] - Log[-(1/(I a + b))] + Log[1/(I a + b)])/
E^(I a + b)^2)))
Although this involves I, numerical evaluation for real a and b (to any
desired precision, followed by Chop, yields the same answer as
nint[a_, b_, opts___] := NIntegrate[Exp[-t^2]/(a^2 + (b - t)^2),
{t, -Infinity, Infinity}, opts]
For example,
nint[2,3]
0.14563
int[2.,3.]//Chop
0.14563
nint[2,3,WorkingPrecision->30]
0.14562973135699681308
int[2`30,3`30]//Chop
0.1456297313569968130774404707
Cheers,
Paul
_______________________________________________________________________
Paul Abbott Phone: 61 8 6488 2734
School of Physics, M013 Fax: +61 8 6488 1014
The University of Western Australia (CRICOS Provider No 00126G)
AUSTRALIA http://physics.uwa.edu.au/~paul