how to compute roots of a cubic function with minimal truncationerrors?

Discussion in 'C++' started by Peng Yu, Sep 10, 2008.

  1. Peng Yu

    Peng Yu Guest

    Hi,

    I have the following program which computes roots of a cubic function.
    The solution is sensitive to the type, which is due to the truncation
    error. 'long double T' gives three solutions, and 'typedef double T'
    gives one solutions. The correct number of solutions should be two, 1
    and 2.

    I know there is some trick to reduce the chance of under or overflow.
    For example, std::abs(z) shall be implemented as
    |x| * sqrt(1 + (y/x)*(y/x)) if |x|>|y|
    and
    |y| * sqrt(1 + (x/y)*(x/y)) if |y|>|x|,
    where z is a complex number, and x and y are its real and complex
    parts.

    I'm wondering what trick can be played to reduce the truncation error
    when solving the cubic polynomial equations.

    BTW, I use the algorithm is shown at http://www.hawaii.edu/suremath/jrootsCubic.html

    Thanks,
    Peng

    #include <vector>
    #include <complex>
    #include <cmath>
    #include <iostream>

    template <typename T>
    std::vector<T> roots_of_cubic_function(const T &a2, const T &a1, const
    T &a0) {

    const T one = static_cast<T>(1);
    const T three = static_cast<T>(3);

    T q = one / 3 * a1 - one / 9 * a2 * a2;
    T r = one / 6 * (a1 * a2 - three * a0) - one / 27 * std::pow(a2, 3);

    T Delta = std::pow(q, 3) + r * r;
    std::cout << "Delta = " << Delta << std::endl;

    std::vector<T> v;
    if(Delta >= T()) {
    T s1 = std::pow(r + sqrt(Delta), one/3);
    T s2 = std::pow(r - sqrt(Delta), one/3);
    v.push_back(s1 + s2 - a2 / 3);
    if(Delta == T()) {
    v.push_back(-.5 * (s1 + s2) - a2 / 3);
    }
    }
    else {
    std::complex<T> temp = sqrt(std::complex<T>(Delta));
    std::complex<T> s1 = std::pow(r + temp, one/3);
    std::complex<T> s2 = std::pow(r - temp, one/3);
    const T minus_half = - static_cast<T>(1)/2;
    v.push_back((s1 + s2 - a2 / 3).real());
    v.push_back((minus_half * (s1 + s2) - a2 / 3 + std::complex<T>(0,
    sqrt(three)/2) * (s1 - s2)).real());
    v.push_back((minus_half * (s1 + s2) - a2 / 3 - std::complex<T>(0,
    sqrt(three)/2) * (s1 - s2)).real());
    }
    return v;
    }

    int main () {
    //typedef long double T;
    typedef double T;
    const T a2 = -4;
    const T a1 = 5;
    const T a0 = -2;

    std::vector<T> v = roots_of_cubic_function<T>(a2, a1, a0);

    std::cout << "Solutions:" << std::endl;
    for(std::vector<T>::const_iterator it = v.begin(); it != v.end(); ++
    it) {
    T x = *it;
    T f = ((x + a2) * x + a1) * x + a0;
    std::cout << x << " " << f << std::endl;
    }
    }
     
    Peng Yu, Sep 10, 2008
    #1
    1. Advertising

  2. Re: how to compute roots of a cubic function with minimal truncation errors?

    [Hi Peng, followup to topical group sci.math.num-analysis]

    From: Peng Yu <>
    Date: Tue, 9 Sep 2008 21:24:10 -0700 (PDT)

    Hi,

    I have the following program which computes roots of a cubic
    function.
    The solution is sensitive to the type, which is due to the truncation
    error. 'long double T' gives three solutions, and 'typedef double T'
    gives one solutions. The correct number of solutions should be two, 1
    and 2.

    I know there is some trick to reduce the chance of under or overflow.
    For example, std::abs(z) shall be implemented as
    |x| * sqrt(1 + (y/x)*(y/x)) if |x|>|y|
    and
    |y| * sqrt(1 + (x/y)*(x/y)) if |y|>|x|,
    where z is a complex number, and x and y are its real and complex
    parts.

    I'm wondering what trick can be played to reduce the truncation error
    when solving the cubic polynomial equations.

    BTW, I use the algorithm is shown at
    http://www.hawaii.edu/suremath/jrootsCubic.html

    Thanks,
    Peng

    #include <vector>
    #include <complex>
    #include <cmath>
    #include <iostream>

    template <typename T>
    std::vector<T> roots_of_cubic_function(const T &a2, const T &a1,
    const
    T &a0) {

    const T one = static_cast<T>(1);
    const T three = static_cast<T>(3);

    T q = one / 3 * a1 - one / 9 * a2 * a2;
    T r = one / 6 * (a1 * a2 - three * a0) - one / 27 * std::pow(a2,
    3);

    T Delta = std::pow(q, 3) + r * r;
    std::cout << "Delta = " << Delta << std::endl;

    std::vector<T> v;
    if(Delta >= T()) {
    T s1 = std::pow(r + sqrt(Delta), one/3);
    T s2 = std::pow(r - sqrt(Delta), one/3);
    v.push_back(s1 + s2 - a2 / 3);
    if(Delta == T()) {
    v.push_back(-.5 * (s1 + s2) - a2 / 3);
    }
    }
    else {
    std::complex<T> temp = sqrt(std::complex<T>(Delta));
    std::complex<T> s1 = std::pow(r + temp, one/3);
    std::complex<T> s2 = std::pow(r - temp, one/3);
    const T minus_half = - static_cast<T>(1)/2;
    v.push_back((s1 + s2 - a2 / 3).real());
    v.push_back((minus_half * (s1 + s2) - a2 / 3 + std::complex<T>(0,
    sqrt(three)/2) * (s1 - s2)).real());
    v.push_back((minus_half * (s1 + s2) - a2 / 3 - std::complex<T>(0,
    sqrt(three)/2) * (s1 - s2)).real());
    }
    return v;
    }

    int main () {
    //typedef long double T;
    typedef double T;
    const T a2 = -4;
    const T a1 = 5;
    const T a0 = -2;

    std::vector<T> v = roots_of_cubic_function<T>(a2, a1, a0);

    std::cout << "Solutions:" << std::endl;
    for(std::vector<T>::const_iterator it = v.begin(); it != v.end();
    ++
    it) {
    T x = *it;
    T f = ((x + a2) * x + a1) * x + a0;
    std::cout << x << " " << f << std::endl;
    }
    }


    --
    Quidquid latine scriptum est, altum videtur.
     
    Martin Eisenberg, Sep 10, 2008
    #2
    1. Advertising

  3. Peng Yu

    Raymond Toy Guest

    Re: how to compute roots of a cubic function with minimal truncation errors?

    >>>>> "Peng" == Peng Yu <> writes:

    Peng> Hi,
    Peng> I have the following program which computes roots of a cubic function.
    Peng> The solution is sensitive to the type, which is due to the truncation
    Peng> error. 'long double T' gives three solutions, and 'typedef double T'
    Peng> gives one solutions. The correct number of solutions should be two, 1
    Peng> and 2.

    Peng> I know there is some trick to reduce the chance of under or overflow.
    Peng> For example, std::abs(z) shall be implemented as
    Peng> |x| * sqrt(1 + (y/x)*(y/x)) if |x|>|y|
    Peng> and
    Peng> |y| * sqrt(1 + (x/y)*(x/y)) if |y|>|x|,
    Peng> where z is a complex number, and x and y are its real and complex
    Peng> parts.

    Peng> I'm wondering what trick can be played to reduce the truncation error
    Peng> when solving the cubic polynomial equations.

    Peng> BTW, I use the algorithm is shown at http://www.hawaii.edu/suremath/jrootsCubic.html

    Look at the formula for s1 and s2. You might want to factor the r out
    from sqrt(q^3+r^2) to abs(r)*sqrt(1+q^3/r^2). This might give you
    better accuracy for s1 and s2.

    Or choose one of the real roots, use Newton iteration to refine it,
    and then solve the resulting quadratic separately. This requires a
    bit of care too.

    Ray
     
    Raymond Toy, Sep 12, 2008
    #3
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. JKop

    Cubic Root

    JKop, Apr 24, 2004, in forum: C++
    Replies:
    8
    Views:
    6,218
    Gunnar G
    Apr 29, 2004
  2. Will Ware
    Replies:
    0
    Views:
    444
    Will Ware
    Nov 21, 2005
  3. jitasi
    Replies:
    1
    Views:
    753
    Terry Reedy
    Mar 4, 2007
  4. Peng Yu
    Replies:
    4
    Views:
    515
    Michael DOUBEZ
    Sep 10, 2008
  5. PerlFAQ Server
    Replies:
    0
    Views:
    269
    PerlFAQ Server
    Feb 2, 2011
Loading...

Share This Page