Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > C++ > how to compute roots of a cubic function with minimal truncationerrors?

Reply
Thread Tools

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

 
 
Peng Yu
Guest
Posts: n/a
 
      09-10-2008
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:ow(a2, 3);

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

std::vector<T> v;
if(Delta >= T()) {
T s1 = std:ow(r + sqrt(Delta), one/3);
T s2 = std:ow(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:ow(r + temp, one/3);
std::complex<T> s2 = std:ow(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;
}
}
 
Reply With Quote
 
 
 
 
Martin Eisenberg
Guest
Posts: n/a
 
      09-10-2008
[Hi Peng, followup to topical group sci.math.num-analysis]

From: Peng Yu <(E-Mail Removed)>
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:ow(a2,
3);

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

std::vector<T> v;
if(Delta >= T()) {
T s1 = std:ow(r + sqrt(Delta), one/3);
T s2 = std:ow(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:ow(r + temp, one/3);
std::complex<T> s2 = std:ow(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.
 
Reply With Quote
 
 
 
 
Raymond Toy
Guest
Posts: n/a
 
      09-12-2008
>>>>> "Peng" == Peng Yu <(E-Mail Removed)> 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
 
Reply With Quote
 
 
 
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
Re: How include a large array? Edward A. Falk C Programming 1 04-04-2013 08:07 PM
FAQ 4.43 How do I compute the difference of two arrays? How do I compute the intersection of two arrays? PerlFAQ Server Perl Misc 0 02-02-2011 05:00 AM
are there any "Piecewise Cubic Hermite Interpolating Polynomial" in python jitasi Python 1 03-04-2007 07:20 AM
Piecewise-cubic lookup table generator Will Ware Python 0 11-21-2005 10:19 PM
Cubic Root JKop C++ 8 04-29-2004 08:32 PM



Advertisments