Homec4science
Diffusion lammps (lammm-devel) Browse

lammps/lib/atc/PolynomialSolver.cpp
8d47d9dc3965lammm-devel

Last Change
Raw File

PolynomialSolver.cpp
View Options

#include "PolynomialSolver.h"
#include <limits>
#include <cmath>
#include <iostream>
#include "ATC_Error.h"
namespace ATC {
// Utility functions used by solvers, but not globally accessible.
static const double PI_OVER_3 = acos(-1.0)*(1.0/3.0);
static bool is_zero(double x)
{
static double GT_ZERO = 1.0e2*std::numeric_limits<double>::epsilon();
static double LT_ZERO = -GT_ZERO;
return x>LT_ZERO && x<GT_ZERO;
}
static double sign(double x)
{
static double s[] = {-1.0,1.0};
return s[x>0];
}
// Linear solver
int solve_linear(double c[2], double x0[1])
{
if (c[1] == 0) return 0; // constant function
*x0 = -c[0] / c[1];
return 1;
}
// Quadratic solver
int solve_quadratic(double c[3], double x0[2])
{
if (is_zero(c[2])) return solve_linear(c, x0);
const double ainv = 1.0/c[2]; // ax^2 + bx + c = 0
const double p = 0.5 * c[1] * ainv; // -b/2a
const double q = c[0] * ainv; // c/a
double D = p*p-q;
if (is_zero(D)) { // quadratic has one repeated root
x0[0] = -p;
return 1;
}
if (D > 0) { // quadratic has two real roots
D = sqrt(D);
x0[0] = D - p;
x0[1] = -D - p;
return 2;
}
return 0; // quadratic has no real roots
}
// Cubic solver
int solve_cubic(double c[4], double x0[3])
{
int num_roots;
if (is_zero(c[3])) return solve_quadratic(c, x0);
// normalize to x^3 + Ax^2 + Bx + C = 0
const double c3inv = 1.0/c[3];
const double A = c[2] * c3inv;
const double B = c[1] * c3inv;
const double C = c[0] * c3inv;
// substitute x = t - A/3 so t^3 + pt + q = 0
const double A2 = A*A;
const double p = (1.0/3.0)*((-1.0/3.0)*A2 + B);
const double q = 0.5*((2.0/27.0)*A*A2 - (1.0/3.0)*A*B + C);
// Cardano's fomula
const double p3 = p*p*p;
const double D = q*q + p3;
if (is_zero(D)) {
if (is_zero(q)) { // one triple soln
x0[0] = 0.0;
num_roots = 1;
}
else { // one single and one double soln
const double u = pow(fabs(q), 1.0/3.0)*sign(q);
x0[0] = -2.0*u;
x0[1] = u;
num_roots = 2;
}
}
else {
if (D < 0.0) { // three real roots
const double phi = 1.0/3.0 * acos(-q/sqrt(-p3));
const double t = 2.0 * sqrt(-p);
x0[0] = t * cos(phi);
x0[1] = -t * cos(phi + PI_OVER_3);
x0[2] = -t * cos(phi - PI_OVER_3);
num_roots = 3;
}
else { // one real root
const double sqrt_D = sqrt(D);
const double u = pow(sqrt_D + fabs(q), 1.0/3.0);
if (q > 0) x0[0] = -u + p / u;
else x0[0] = u - p / u;
num_roots = 1;
}
}
double sub = (1.0/3.0)*A;
for (int i=0; i<num_roots; i++) x0[i] -= sub;
return num_roots;
}
// solve ode with polynomial source : y'n + a_n-1 y'n-1 + ... = b_n x^n +...
void integrate_ode(double x,
int na, double * a, double * y0, double * y, int nb, double *b )
{
if (na == 2) {
// particular
if ( a[1] == 0) {
if ( a[0] == 0) {
y[0] = y0[0]+y0[1]*x;
y[1] = y0[1];
}
else {
double c = sqrt(a[0]);
y[0] = y0[0]*cos(c*x)+y0[1]/c*sin(c*x);
y[1] = -c*y0[0]*cos(c*x)+y0[1] *sin(c*x);
}
}
else {
// use solve_quadratic
throw ATC_Error("not yet supported");
}
// homogenous
double c = 1.;
double z = x;
int j = 2;
for (int i = 0; i < nb; i++,j++) {
y[1] += j*c*z;
c /= j;
z *= x;
y[0] += c*z;
}
}
else throw ATC_Error("can only integrate 2nd order ODEs currently");
}
}
c4science · Help