What are good methods for detecting parabolic components and Siegel disk components in the Fatou set of a rational function? Given a rational function $f$ on the Riemann sphere, I would like to answer the question: Does the Fatou set of the function, $F(f)$, contain any parabolic components or Siegel disk components?  Additionally, if it turns out that $F(f)$ contains parabolic components or Siegel disk components I would like to visualize those components, perhaps by highlighting their boundaries in a plot of the Julia set of $f$.  My question is:

Is there a known computational method for reliably detecting and describing the parabolic components and Siegel disk components of the Fatou set of an arbitrary rational function?

I am aware that a Siegel disk component will contain an indifferent periodic point of $f$ and that a parabolic component will have an indifferent periodic point of $f$ on its boundary, so a preliminary step in the method I seek might be to merely identify the indifferent periodic points of $f$. I have thought of a method for detecting these points, but I'm sure it is unreliable.  The method is to take a large number of sample points, $z_0$, and use any cycle detection algorithm (say tortoise and hare) on the sequence $f^n(z_0)$.  Then if that process detects a cycle of period $m$ at $w$ (equal to one of $z_0$, $f(z_0)$, $f^2(z_0)$, etc.), to numerically measure the magnitude of $(f^m)'(w)$ to see if that magnitude is tolerably close to 1.  My secondary question is:

Is there  a known efficient and reliable computational method for identifying the indifferent periodic points of an arbitrary rational function?

Update: I have made some progress by following Alexandre Eremenko's suggestion to examine the orbits of the critical points.  I would like to share an example application of this technique for identifying (rational) neutral cycles.  
Consider the family of functions $f:z\to z - p(z)/p'(z)$, describing Newton-Raphson iteration on $p$, where $p$ is a cubic polynomial with simple roots.  An element $f$ of this family has three attracting fixed points in $\mathbb{C}$, the roots of $p$, and is conjugate by a Möbius transformation to a function with an attracting fixed point at -1, an attracting fixed point at 1, and an attracting fixed point at $\lambda\in\mathbb{C}$.  The following image is a visualization of types of cycles occurring in $f_\lambda$, for $\lambda$ in the axis aligned square of side length 0.04 with top left corner $-0.834 + 0.861 i$.  The only candidate for a critical point of $f_\lambda$ is $\alpha = \lambda/3$.

In the image a positive green color component at a point $\lambda$ indicates that $f_\lambda^n(\alpha)$ converges to one of the fixed points (1, -1, or $\lambda$).  Brighter green indicates slower convergence.  The green regions are typical of Newton basin renderings.  Points $\lambda$ having no green color component, but positive red and blue color components, indicate that $f_\lambda^n(\alpha)$ converges to a cycle of period greater than 1.  In this case the cycle may be neutral or attracting.  Brighter purple indicates slower convergence, so the bright purple regions indicate values of $\lambda$ such that $f_\lambda$ probably has neutral periodic points and $F(f_\lambda)$ might have parabolic components or Siegel disk components.
 A: It depends on how your function is given, and what exactly you mean by "detecting".
Anyway, you don't have to look at the "large number of sample points $z_0$.
It is enough to look at the critical points. In principle, all neutral rational cycles
and Siegel discs can be found from the behavior of the orbits of critical points.
Every attracting cycle and every 
rational neutral cycle attracts some critical point (Fatou's theorem).
Find many points of the orbit of every critical point. If you see an orbit converging to a cycle,
you can easily tell an attracting cycle from a neutral rational one: convergence to an attracting
cycle is geometric, while convergence to a neutral rational cycle is EXTREMELLY slow.
Detecting a Siegel disc, it is much harder. It is known that the orbits of critical points
are dense on the boundary of a Siegel disc. But this is of little help.
But at least a picture of the orbits of critical points can give you the idea of what the preiod of your Siegel disc could
possibly be. Then you just find the fixed points of the corresponding iterate and see whether
any of them is neutral.
Actually, if your function is given numerically, it is unlikely that it has a neutral point,
so there is nothing to detect...
A: For detecting Siegel disc around fixed point in a plot of the Julia set of f one can use average discrete velocity of orbit.

(source)
It is described here
https://commons.wikimedia.org/wiki/File%3AGolden_Mean_Quadratic_Siegel_Disc_Speed.png
I have learned this method from Chris King
http://www.dhushara.com/DarkHeart/DarkHeart.htm
To highlight boundaries  of Siegel disk components in a plot of the Julia set of $f$ one can use this method : 

(source)
It works for cases when critical point is on the boundary of Siegel disc component.
See :
Building blocks for Quadratic Julia sets by : Joachim Grispolakis, John C. Mayer and Lex G. Oversteegen Journal: Trans. Amer. Math. Soc. 351 (1999), 1171-1201
HTH
Adam
A: "Detecting" whether a neutral cycle exists or not will be difficult in general. 
However, when you do know that there is a parabolic cycle (e.g. because you have constructed your map that way), Braverman shows, as quoted by Adam, that you can compute the Julia set in polynomial time. What should be noted here is that the program used (and the time bound) will depend on the parameter in question. 
For irrationally indifferent fixed points, things can get even more tricky. But when the rotation number is nice (e.g. some quadratic irrational such as the 'golden mean'), you will have a critical point on the boundary of the Siegel disk. So you can draw the boundary of the Siegel disk by plotting the orbit of the critical point.
In parameter spaces, rather than trying to 'detect' neutral cycles, you may wish to draw the boundaries of hyperbolic components using Newton's method. That is, take a point in the hyperbolic component that you are interested in (where there is an attracting cycle), and then find a curve along which the modulus of the multiplier tends to one. Then you will have found an indifferent parameter. Now you can similarly change the argument of the multiplier, again using Newton's method, and trace this curve. Some care is required near "cusps". 
For an example of such a picture, in the family of exponential maps, see Figure 1 in my article with Dierk Schleicher, "Bifurcations in the space of exponential maps" (http://arxiv.org/abs/math/0311480).
A: To visualize parabolic  components by highlighting their boundaries in a plot of the Julia set of f one can use :


*

*giant steps method by Mark Braverman ( http://arxiv.org/abs/math.DS/0505036 )

*modified distance etimation method (image, src code and comments :
http://commons.wikimedia.org/wiki/File:Parabolic_Julia_set_for_internal_angle_1_over_5.png)

*combined methods : http://commons.wikimedia.org/wiki/File:Parabolic_Julia_set_for_internal_angle_1_over_4.png
Some numerical problems ( lazy dynamics) are described here :
http://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/parabolic
HTH
