comp.lang.idl-pvwave archive
Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado

Home » Public Forums » archive » Cubic root finding on a grid
Show: Today's Messages :: Show Polls :: Message Navigator
E-mail to friend 
Switch to threaded view of this topic Create a new topic Submit Reply
Cubic root finding on a grid [message #68737] Fri, 20 November 2009 08:03 Go to next message
Luds is currently offline  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 Go to previous message
Lajos Foldy is currently offline  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 Go to previous message
wlandsman is currently offline  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 Go to previous message
Lajos Foldy is currently offline  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
  Switch to threaded view of this topic Create a new topic Submit Reply
Previous Topic: Re: Finding distance with longitude and latitude
Next Topic: All logical unit currently in use

-=] Back to Top [=-
[ Syndicate this forum (XML) ] [ RSS ] [ PDF ]

Current Time: Wed Oct 08 11:51:27 PDT 2025

Total time taken to generate the page: 0.00598 seconds