Cubic root finding on a grid [message #68737] |
Fri, 20 November 2009 08:03  |
Luds
Messages: 7 Registered: September 2009
|
Junior Member |
|
|
Hi IDL'ers,
I'm trying to find an efficient method for estimating the eigenvalues
of a tensor which is defined on a cubic (NxNxN) grid. The tensor is a
simple 3x3 (symmetric) matrix defined at each of the N^3 grid points,
so the eigenvalues at a given grid point can be calculated by (for
example) the IDL routine IMSL_ZEROPOLY simply by finding the roots of
the cubic polynomial defined by the matrix at that point - but this
routine doesn't handle array's of coefficients, so it has to be
evaluated on the grid point by point.
My grids are up to 512^3, so using a for-loop to compute the
eigenvalues at each node is rather slow. Does anyone know of any
adaptations of the IMSL_ZEROPOLY routine that can work on a grid of
3X3 matricies (or on an array of polynomial coefficients)?
Any suggestions would help.
Best,
Ludlow
|
|
|
Re: cubic root [message #83959 is a reply to message #68737] |
Tue, 16 April 2013 09:44  |
Lajos Foldy
Messages: 176 Registered: December 2011
|
Senior Member |
|
|
On Tuesday, April 16, 2013 5:19:00 PM UTC+2, fawltyl...@gmail.com wrote:
> (-1)^1.0/3 gives -0.333333 for me.
>
> You have to specify that you want complex roots:
>
> IDL> help, complex(-1.0)^(1.0/3)
> <Expression> COMPLEX = ( 0.500000, 0.866025)
>
Also, here is a function for the real cubic root:
function cubic_root, x
return, x ge 0 ? x^(1/3.0) : -(-x)^(1/3.0)
end
(Change 3.0 to 3.0d for double precision x.)
regards,
Lajos
|
|
|
Re: cubic root [message #83963 is a reply to message #68737] |
Tue, 16 April 2013 08:27  |
wlandsman
Messages: 743 Registered: June 2000
|
Senior Member |
|
|
You can also use FZ_ROOTS to get all 3 roots to the equation x^3 + 1 = 0
IDL> print,fz_roots([1,0d,0,1.0d])
( -1.0000000, 0.0000000)( 0.50000000, -0.86602540)( 0.50000000, 0.86602540)
or use http://idlastro.gsfc.nasa.gov/ftp/contrib/freudenreich/cuber oot.pro
On Tuesday, April 16, 2013 11:08:12 AM UTC-4, fd_...@mail.com wrote:
> Hi all
>
>
>
> I have a question about cubic roots. The cubic roots of a negative number there exists. The cubir root of -1 is -1, 0f -2 is -1.25 and so on why when I typed in IDL to calculate the cubic root of a negative value it print out NaN which means Not a Number?
>
>
>
> E.g. I typed: (-1)^1.0/3 and it prints out NaN
>
>
>
> Is there any way to fixed it and don't get NaNs but get the actual value?
>
>
>
> Best Wishes
>
> Maria
|
|
|
Re: cubic root [message #83964 is a reply to message #68737] |
Tue, 16 April 2013 08:19  |
Lajos Foldy
Messages: 176 Registered: December 2011
|
Senior Member |
|
|
On Tuesday, April 16, 2013 5:08:12 PM UTC+2, fd_...@mail.com wrote:
> Hi all
>
>
>
> I have a question about cubic roots. The cubic roots of a negative number there exists. The cubir root of -1 is -1, 0f -2 is -1.25 and so on why when I typed in IDL to calculate the cubic root of a negative value it print out NaN which means Not a Number?
>
>
>
> E.g. I typed: (-1)^1.0/3 and it prints out NaN
>
>
>
> Is there any way to fixed it and don't get NaNs but get the actual value?
>
>
>
> Best Wishes
>
> Maria
(-1)^1.0/3 gives -0.333333 for me.
You have to specify that you want complex roots:
IDL> help, complex(-1.0)^(1.0/3)
<Expression> COMPLEX = ( 0.500000, 0.866025)
(and note that ^ has higher precedence than /).
regards,
Lajos
|
|
|