Dear all,
I have a question. I am a new Root User and I have to implement a code where I need to calculate complex function on complex number. How can I define a complex function in Root? I have tried with TComplex with something like
TF1 *fp = new TF1 (“fp”, “cos(x)+x**TComplex::I()”, 0., 10.);
or
fa1=new TF1(“fa1”, “TComplex(cos(x),x)”, 0, 10);
but it didn’t work. There are rapid method or I need to write a new code?
Thank you,
Eleonora
[code]#include
#include
#include
std::complex MyFunc(double x) // cos(x) + x^i (note: x >= 0)
{
std::complex v(0.0);
if (x == 0.0) return std::complex(1.0, 0.0);
if (x < 0.0) std::cout << “Error : ‘x’ MUST NOT be negative!” << std::endl;
else v = std::complex( std::cos(x) + std::cos(std::log(x)),
std::sin(std::log(x)) );
return v;
}[/code]
Hi,
Thank you very much for the reply. For know if I have understand, in this way you have used a mathematical identity
cos(x)+ix=cos(x)+cos(log(x))+isin(log(x))= cos(x)+exp(i^log(x))?
and so I can use only x>=0 because of the log, right?
In your original post, in the “fp” function, I can see “xTComplex::I()", that is "xi”, that is “x^i”.
Hi,
Sorry, in copy and paste I have posted the wrong function!
Thank you for the reply!
Eleonora
[code]#include
#include
std::complex MyFunc(double x) // cos(x) + x*i
{ return std::complex( std::cos(x), x ); }[/code]
Thank you again!
Probably a stupid question. What if I want the same function but on complex number? I have tried
std::complex MyFunc3(complex x)
{ return std::complex( std::cos(x), x ); }
but it doesn’t work: “Symbol null is not defined in current scope”…
[code]#include
std::complex MyFunc(std::complex x) // cos(x) + x*i
{ return std::cos(x) + x * std::complex(0.0, 1.0); }[/code]
See also:
- Complex numbers library - header reference @ cplusplus.com
- C numerics library - header reference @ cplusplus.com
Thank you =)
I need to definite the Hankel function, that can be written with the Bessel function JBessel and the Neuman function yNeuman:
h1_n(x)=j_n(x)+i*y_n(x),
so h1 is my complex function.
I have written code in different way but I can’t find the right mode.
In this case I have tried only with the first component,
std::complex Hankel0(std::complex x)
{
gSystem->Load(“libMathMore”);
TF1* jBessel0= new TF1(“J_0”, “ROOT::Math::sph_bessel(0,x)”, 0, 10);
TF1* yBessel0= new TF1(“J_0”, “ROOT::Math::sph_neumann(0,x)”, 0, 10);
return jBessel0 + yBessel0* std::complex(0.0, 1.0);
}
but it didn’t work. Where I am wrong? Thank you for the patience
Create a “Hankel0.cxx” file: [code]#include “Math/SpecFunc.h”
#include
#include
// sph_bessel(0, x) + i * sph_neumann(0, x)
std::complex Hankel0(double x) // note: x > 0
{
std::complex v(0.0);
if (x <= 0.0) std::cout << “Error : ‘x’ MUST be positive!” << std::endl;
else v = std::complex( ROOT::Math::sph_bessel(0, x),
ROOT::Math::sph_neumann(0, x) );
return v;
}[/code] and then try:
root [0] #include
root [1] gSystem->Load(“libMathMore”);
root [2] .L Hankel0.cxx++
root [3] std::cout << Hankel0(1.0) << std::endl;
Thank =)
Hi,
I have defined the PsiBessel jBessel*x function that works on real numbers
double psiBessel0(double x)
{return ROOT::Math::sph_bessel(0, x)*x;}
but I need also a PsiBessel that works on complex number, so I defined
std::complex psiBessel0C(std::complex x)
{return ROOT::Math::sph_bessel(0, x)* x * std::complex(0.0, 1.0);}
but It didn’t work and I don’t understand why.
And when you check the description of the “ROOT::Math::sph_bessel”, can you find any function prototype that accepts “complex” variables (and/or returns “complex” results)?
I find this
spherical Bessel functions of the first kind
double sph_bessel(unsigned n, double x)
{
return gsl_sf_bessel_jl(n, x);
}
So you think that maybe this function in root is defined only on real numbers?
See: The GNU Scientific Library - Reference Manual - 7.5.5 Regular Spherical Bessel Functions
Note that the “Bessel Function of the First Kind” is NOT the same as the “Spherical Bessel Function of the First Kind”.
Similarly, the “Bessel Function of the Second Kind (a.k.a. Neumann function or Weber function)” is NOT the same as the “Spherical Bessel Function of the Second Kind”.
Then, the “Hankel Function” and/or the “Hankel Function of the First Kind (a.k.a. Bessel function of the third kind or Weber function)” is NOT the same as the “Spherical Hankel Function of the First Kind”.
(Note also that, in the web links given above, the “spherical” functions of the orders 0 … 3 are explicitly expressed in terms of simple elementary functions, even for “complex” arguments.)
Thank for the link: I have checked, I need the Spherical Bessel of the first and second kind, so I can write also the Spherical Hankel of the first kind. Sorry for my ignorance but I known very little of root or C++: so these in ROOT::Math are defined on double so I can use them only on real number., and If I need to use them on complex I have to define the function (for example using the spherical definition of orders 0 … 3 expressed in elementary functions), Right?
I would try to proceed this way.
One needs to take care of some “pathological” cases when Re(x) and/or Im(x) are equal to 0 (for example one needs to identify terms in form “sin(some_real_value)/some_real_value” and “1/some_real_value”, which need to be taken care of -> you could also rewrite all equations using explicit separation of their “real” and “imaginary” parts, perform calculations on two “real values” and then “combine” the results into a single “complex” number).
Maybe Lorenzo has any better idea.
[code]#include
#include
#include
double SinXOverX(const double x) // = sin(x) / x
{
if (std::fabs(x) < 1e-7) return 1.0 - x * x / 6.0;
return std::sin(x) / x;
}
double CosXOverX(const double x) // = cos(x) / x
{
if (std::fabs(x) < 1e-7) {
if (std::fabs(x) < 1e-300) {
std::cout << “CosXOverX : calculations went amok?” << std::endl;
if (x < 0.0) return -1e300;
return 1e300;
}
return 1.0 / x - x / 2.0;
}
return std::cos(x) / x;
}
// Spherical Bessel Function of the First Kind (order 0)
// http://mathworld.wolfram.com/SphericalBesselFunctionoftheFirstKind.html
std::complex SBesselF0(std::complex x) // = sin(x) / x
{
if (x.imag() == 0.0) return std::complex(SinXOverX(x.real()), 0.0);
return std::sin(x) / x;
}
// PsiSBesselF0 = x * SBesselF0(x)
std::complex PsiSBesselF0(std::complex x) // = sin(x)
{
return std::sin(x);
}
// Spherical Bessel Function of the Second Kind (order 0)
// http://mathworld.wolfram.com/SphericalBesselFunctionoftheSecondKind.html
std::complex SBesselS0(std::complex x) // = - cos(x) / x
{
if (x.imag() == 0.0) return std::complex(-CosXOverX(x.real()), 0.0);
return (-std::cos(x) / x);
}
// PsiSBesselS0 = x * SBesselS0(x)
std::complex PsiSBesselS0(std::complex x) // = - cos(x)
{
return (-std::cos(x));
}
// Spherical Hankel Function of the First Kind (order 0)
// http://mathworld.wolfram.com/SphericalHankelFunctionoftheFirstKind.html
std::complex SHankelF0(std::complex x) // = -i * exp(i * x) / x
{
static const std::complex I(0.0, 1.0);
if (x.imag() == 0.0) return std::complex( SinXOverX(x.real()),
-CosXOverX(x.real()) );
return std::conj(I) * std::exp(I * x) / x;
}
// PsiSHankelF0 = x * SHankelF0(x)
std::complex PsiSHankelF0(std::complex x) // = -i * exp(i * x)
{
static const std::complex I(0.0, 1.0);
return std::conj(I) * std::exp(I * x);
}[/code]