By using this site, you agree to our updated Privacy Policy and our Terms of Use. Manage your Cookies Settings.
448,563 Members | 1,205 Online
Bytes IT Community
+ Ask a Question
Need help? Post your question and get tips & solutions from a community of 448,563 IT Pros & Developers. It's quick & easy.

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

P: n/a
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<Troots_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<Tv;
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<Ttemp = sqrt(std::complex<T>(Delta));
std::complex<Ts1 = std::pow(r + temp, one/3);
std::complex<Ts2 = 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<Tv = 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;
}
}
Sep 10 '08 #1
Share this Question
Share on Google+
2 Replies


P: n/a
[Hi Peng, followup to topical group sci.math.num-analysis]

From: Peng Yu <Pe*******@gmail.com>
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<Troots_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<Tv;
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<Ttemp = sqrt(std::complex<T>(Delta));
std::complex<Ts1 = std::pow(r + temp, one/3);
std::complex<Ts2 = 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<Tv = 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.
Sep 10 '08 #2

P: n/a
>>>>"Peng" == Peng Yu <Pe*******@gmail.comwrites:

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

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

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

PengBTW, 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
Sep 12 '08 #3

This discussion thread is closed

Replies have been disabled for this discussion.