Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.


Regula falsi vs. secant method

  • Mihai

    Hey Rudiger,

    I didn't receive an email, I didn't know you answered to this question.
    Now that I see it, I respond.

    I understand the idea, I still look at the code (I'm confused about the virtual F method in your code, that must be defined in the case of every range (if I'm not mistaken).

    As I see it, in those cases of 'higher-order contacts' as you named them, the planet speeds of the sides of the interval are opposites, right?
    When this happens, it may worth it looking, in little increments, if the tangent (first derivative?) has any chance to reach desired longitude.
    When this is impossible, then there's no zero there.

    This was the special case which took me down whenever I wanted to implement this algorithm myself.
    but you say there are small chances (I see them small too, but, it was enough to take me down)


  • Mihai

    Actually, when trying to implement an algorithm to search simple patterns, I didn't followed the 'traditional way' of thinking, with secant method.

    I tried something (inovative?) trying to use average speed of the planet, speed sign (less or over 0), and distance sign to try to jump in the future, then to analyze where I am, then to jump again forward or backward until I'm in the accepted error near desired longitude.
    Now, when I write here, I think other fixed parameters would be useful too (like maximum and minimum values for each parameter for each planet)

    I simply don't like this method of daily jumps in a range.

    But it's the only one it works, so far.

    In your code it works quite fast. Maybe you did some improvements with incr variable.
    I saw that swiss ephemeris can do something like 10000 calculations/second. So maybe this method of jumping in increments is ok after all.

    And I have something to say about searching patterns and daily jumps (or incr jumps):
    What I say here is about my method of "adaptive jumping"
    As I imagined, when trying to find patterns, I saw that first I should find all cases only with the biggest planet first. Then continue with the next smaller planet, then the next, until I find a pattern. In this way it will be quite fast, instead of making daily jumps (or incr jumps) to search the perfect match.

    I really would love your response about this way of 'adaptive jumping'

    I'm aware there are big wholes in the picture I see, but I remained at the  level of just picturing this process of searching.


  • Dear  Mihai,

    I looked around a bit how this gap could be filled. I found an interesting method in section 2.6.2 of  Gisela Engeln-Müllges / Frank Uhlig "Numerical algorithms with C", volume 1, available completely online in Google books

    It can be proven for a sufficiently wide class of functions that, if f has a higer order zero at x, then this x is a first-order zero of the function h, defined as 

    h(x) = f(x)^2 / (f(x+f(x)) - f(x))

    Thus, the standard regula falsi applied to h instead of f, gives the zero of f(x). This is called the "modified regula falsi". By a wise usage of inheritance, the existing class Valuation could be reused to afford this alternative search.

    Of course, the performance of AstroPatterns would deteriorate quickly if this second method was called too frequently. For long times, a function like the difference of two planet longitudes, will keep its sign, i.e. stay positive or negative. But the values will be too far away from the zero - even if both planets would have their maximum speed (even in opposite directions!), they could not reach zero inside of the interval considered.

    As you might have seen in the code,  there is a skip() method called first, which allows to skip in the future if the Valuation function is yet too far away from the zero. For usage in the implementations of this skip method, I have a class with the maximum speed values in longitude, declination, latitude, distance,  … I have determined these values with a separate program speed.c which once ran over the complete Swiss Ephemeris time span. I then copied the obtained values into the source code of the class PlanetData.cpp

    The skip() method is near to what you are calling "adaptive jumping", and it's already implemented for several ranges, for example for Solar and Lunar Horoscopes. For aspect ranges, the matter is not that simple, but the values for extremal motion can be helpful there, too.

    Kind regards,