PLUTO
Functions | Variables
quartic.c File Reference

Solve quartic and cubic equations-. More...

#include "pluto.h"
Include dependency graph for quartic.c:

Go to the source code of this file.

Functions

void CheckSolutions (double b, double c, double d, double e, double *z)
 
void PrintSolution (double *z)
 
void QuarticPrintCoeffs (double b, double c, double d, double e)
 
int QuadraticSolve (double a, double b, double c, double *x)
 
void SortArray (double *z, int n)
 
int CubicSolve2 (double b, double c, double d, double z[])
 
int QuarticSolve (double b, double c, double d, double e, double *z)
 
int CubicSolve (double b, double c, double d, double z[])
 
int QuarticSolveNew (double b, double c, double d, double e, double *z)
 

Variables

static int debug_print =0
 

Detailed Description

Solve quartic and cubic equations-.

Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Sept 10, 2012

Definition in file quartic.c.

Function Documentation

void CheckSolutions ( double  b,
double  c,
double  d,
double  e,
double *  z 
)

Definition at line 235 of file quartic.c.

239 {
240  int n;
241  double x, f, fp, fm, fmax;
242 
243 /* -- scale to min and max solution */
244 
245  x = 1.0; fp = e + x*(d + x*(c + x*(b + x)));
246  x = -1.0; fm = e + x*(d + x*(c + x*(b + x)));
247  for (n = 0; n < 4; n++){
248  x = z[n];
249  f = e + x*(d + x*(c + x*(b + x)));
250  print ("> CheckSolutions(): f(z[%d]) = %8.3e\n",n,f);
251  }
252 
253 }
n
static int n
Definition: analysis.c:3
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
GRID::x
double * x
Definition: structs.h:80
menu.c
tuple c
Definition: menu.py:375
fp
FILE * fp
Definition: analysis.c:7

Here is the call graph for this function:

Here is the caller graph for this function:

int CubicSolve ( double  b,
double  c,
double  d,
double  z[] 
)

Solve a cubic equation in the form

\[ z^3 + bz^2 + cz + d = 0 \]

For its purpose, it is assumed that ALL roots are real. This makes things faster.

Parameters
[in]bcoefficient of the cubic
[in]ccoefficient of the cubic
[in]dcoefficient of the cubic
[out]zvector containing the roots of the cubic. Roots should be sorted in increasing order.

Reference:

  • Section 5.6 of Numerical Recipe
Returns
Return 0 on success.

Definition at line 59 of file quartic.c.

80 {
81  double b2;
82  double Q, R;
83  double sQ, arg, theta, cs, sn, p;
84  const double one_3 = 1.0/3.0;
85  const double one_27 = 1.0/27.0;
86  const double sqrt3 = sqrt(3.0);
87 
88  b2 = b*b;
89 
90 /* ----------------------------------------------
91  the expression for f should be
92  Q = c - b*b/3.0; however, to avoid negative
93  round-off making h > 0.0 or g^2/4 - h < 0.0
94  we let c --> c(1- 1.1e-16)
95  ---------------------------------------------- */
96 
97  Q = b2*one_3 - c*(1.0 - 1.e-16); /* = 3*Q, with Q given by Eq. [5.6.10] */
98  R = b*(2.0*b2 - 9.0*c)*one_27 + d; /* = 2*R, with R given by Eq. [5.6.10] */
99 
100 Q = MAX(Q, 0.0);
101 /*
102 if (fabs(Q) < 1.e-18){
103  print ("CubicSolve() very small Q = %12.6e\n",Q);
104  QUIT_PLUTO(1);
105 }
106 if (Q < 0.0){
107  print ("! CubicSolve(): Q = %8.3 < 0 \n",Q);
108  QUIT_PLUTO(1);
109 }
110 */
111 
112 /* -------------------------------------------------------
113  We assume the cubic *always* has 3 real root for
114  which R^2 > Q^3.
115  It follows that Q is always > 0
116  ------------------------------------------------------- */
117 
118  sQ = sqrt(Q)/sqrt3;
119  arg = -1.5*R/(Q*sQ);
120 
121 /* this is to prevent unpleseant situation
122  where both g and i are close to zero */
123 
124  arg = MAX(-1.0, arg);
125  arg = MIN( 1.0, arg);
126 
127 
128  theta = acos(arg)*one_3; /* Eq. [5.6.11], note that pi/3 < theta < 0 */
129 
130  cs = cos(theta); /* > 0 */
131  sn = sqrt3*sin(theta); /* > 0 */
132  p = -b*one_3;
133 
134  z[0] = -sQ*(cs + sn) + p;
135  z[1] = -sQ*(cs - sn) + p;
136  z[2] = 2.0*sQ*cs + p;
137 
138  if(debug_print) {
139  int l;
140  double x, f;
141  print ("===========================================================\n");
142  print ("> Resolvent cubic:\n");
143  print (" g(x) = %18.12e + x*(%18.12e + x*(%18.12e + x))\n", d, c, b);
144  print (" Q = %8.3e\n",Q);
145  print (" arg-1 = %8.3e\n", -1.5*R/(Q*sQ)-1.0);
146 
147  print ("> Cubic roots = %8.3e %8.3e %8.3e\n",z[0],z[1],z[2]);
148  for (l = 0; l < 3; l++){ // check accuracy of solution
149 
150  x = z[l];
151  f = d + x*(c + x*(b + x));
152  print (" verify: g(x[%d]) = %8.3e\n",l,f);
153  }
154 
155  print ("===========================================================\n");
156  }
157 
158  return(0);
159 }
MAX
#define MAX(a, b)
Definition: macros.h:101
MIN
#define MIN(a, b)
Definition: macros.h:104
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
GRID::x
double * x
Definition: structs.h:80
menu.c
tuple c
Definition: menu.py:375
debug_print
static int debug_print
Definition: quartic.c:20

Here is the call graph for this function:

Here is the caller graph for this function:

int CubicSolve2 ( double  b,
double  c,
double  d,
double  z[] 
)

Definition at line 356 of file quartic.c.

357 {
358  double Q, R, R_Q, B, A, sQ, th, tn, sn, cn;
359 
360  Q = (b*b - 3.0*c)/9.0; /* Eq. [5.6.10] */
361  R = (2.0*b*b*b - 9.0*b*c + 27.0*d)/54.0; /* Eq. [5.6.10] */
362  R_Q = R/Q;
363 
364 if (debug_print) print ("R = %8.3e, Q = %8.3e, R^2/Q^3 - 1.0 = %8.3e\n",R,Q,R_Q*R_Q/Q-1.0);
365  if (R_Q*R_Q <= Q){ /*** First case: 3 real roots ***/
366 if (debug_print) print ("> CubicSolve2(): 3 real roots\n");
367  sQ = sqrt(Q);
368  th = acos(R_Q/sQ); /* Eq. [5.6.11] */
369  tn = tan(th/3.0);
370  cn = 1.0/sqrt(1.0 + tn*tn);
371  sn = cn*tn;
372 
373  /* The three roots are given by Eq. (5.6.12) of Numerical Recipe
374  (we switch x3 <-> x2 so that one can verify that x1 < x2 < x3 always):
375 
376  x1 = -2.0*sQ*cos(th/3.0) - aa/3.0;
377  x2 = -2.0*sQ*cos((th - 2.0*CONST_PI)/3.0) - aa/3.0;
378  x3 = -2.0*sQ*cos((th + 2.0*CONST_PI)/3.0) - aa/3.0;
379 
380  To avoid loss of precsion we use
381 
382  cos((th + 2pi/3) = cos(th/3)*cos(2pi/3) - sin(th/3)*sin(2pi/3)
383  cos((th - 2pi/3) = cos(th/3)*cos(2pi/3) + sin(th/3)*sin(2pi/3)
384 
385  Using the fact that the roots must lie in [0,1] and that the spline
386  is monotonically increasing, this is how we pick up the solution: */
387 
388 
389  z[0] = -2.0*sQ*cos(th/3.0) - 1.0/3.0;
390  z[1] = -2.0*sQ*cos((th - 2.0*CONST_PI)/3.0) - b/3.0;
391  z[2] = -2.0*sQ*cos((th + 2.0*CONST_PI)/3.0) - b/3.0;
392 
393  }else{ /*** Second case: one root only ***/
394  print( "! CubicSolve2(): 1 root only, arg = %8.3e !!\n", R_Q*R_Q/Q);
395 
396  A = Q/R;
397  A = fabs(R)*(1.0 + sqrt(1.0 - A*A*Q));
398  A = -DSIGN(R)*pow(A, 1.0/3.0); /* Eq. [5.6.15] */
399  if (A == 0) B = 0.0;
400  else B = Q/A;
401  z[0] = (A + B) - b/3.0; /* Eq. [5.6.17] */
402  z[1] = z[2] = 0.0;
403  }
404 
405 
406  if(debug_print) {
407  int l;
408  double x,f;
409  print ("===========================================================\n");
410  print ("> Resolvent cubic:\n");
411  print (" g(x) = %18.12e + x*(%18.12e + x*(%18.12e + x))\n", d, c, b);
412  print ("> Cubic roots = %8.3e %8.3e %8.3e\n",z[0],z[1],z[2]);
413  for (l = 0; l < 3; l++){ // check accuracy of solution
414 
415  x = z[l];
416  f = d + x*(c + x*(b + x));
417  print (" verify: g(x[%d]) = %8.3e\n",l,f);
418  }
419 
420  print ("===========================================================\n");
421  }
422 
423  return 0;
424 }
DSIGN
#define DSIGN(x)
Definition: macros.h:110
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
GRID::x
double * x
Definition: structs.h:80
menu.c
tuple c
Definition: menu.py:375
debug_print
static int debug_print
Definition: quartic.c:20
CONST_PI
#define CONST_PI
.
Definition: pluto.h:269
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87

Here is the call graph for this function:

void PrintSolution ( double *  z)

Definition at line 346 of file quartic.c.

350 {
351  print ("z = %f %f %f %f\n",z[0],z[1],z[2],z[3]);
352 }
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497

Here is the call graph for this function:

int QuadraticSolve ( double  a,
double  b,
double  c,
double *  x 
)

Solve a quadratic equation in the form

ax^2 + bx + c = 0

Return roots in increasing order

Definition at line 256 of file quartic.c.

265 {
266  double del, sb, q;
267 
268  del = b*b - 4.0*a*c;
269 
270  if (del < 0.0) {
271  print ("! QuadraticSolve(): del = %8.3e, resetting to 0\n",del);
272  del = MAX(del,0.0);
273  }
274 
275  if (del < 0.0) {
276  print ("! QuadraticSolve(): complex roots, del = %8.3e\n",del);
277  return 1; /* No real root */
278  }
279 
280  if (b == 0){
281  x[1] = sqrt(c/a);
282  x[0] = -x[1];
283  }
284 
285  sb = (b > 0.0 ? 1.0:-1.0);
286  q = -0.5*(b + sb*sqrt(del));
287  x[0] = q/a;
288  x[1] = c/q;
289 
290  if (x[0] > x[1]) {
291  double scrh = x[1];
292  x[1] = x[0];
293  x[0] = scrh;
294  }
295 
296  return 0;
297 }
a
static double a
Definition: init.c:135
MAX
#define MAX(a, b)
Definition: macros.h:101
configure.scrh
tuple scrh
Definition: configure.py:200
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
GRID::x
double * x
Definition: structs.h:80
menu.c
tuple c
Definition: menu.py:375
q
static Runtime q
Definition: runtime_setup.c:490

Here is the call graph for this function:

void QuarticPrintCoeffs ( double  b,
double  c,
double  d,
double  e 
)

Definition at line 300 of file quartic.c.

304 {
305  print (" f(x) = %18.12e + x*(%18.12e + x*(%18.12e ",
306  e, d, c);
307  print (" + x*(%18.12e + x*%18.12e)))\n", b, 1.0);
308 
309  print (" b = %18.14e;\n",b);
310  print (" c = %18.14e;\n",c);
311  print (" d = %18.14e;\n",d);
312  print (" e = %18.14e;\n",e);
313 
314  double a2 = -2.0*c;
315  double a1 = c*c + b*d - 4.0*e;
316  double a0 = -(b*c*d - b*b*e - d*d);
317 
318  print (" Resolvent cubic:\n");
319  print (" g(x) = %18.12e + x*(%18.12e + x*(%18.12e + x))\n", a0, a1, a2);
320  double Y = 0.25*b*b + 2.0*d/(b+1.e-40) - c;
321  print (" b^2/4 = %f\n",0.25*b*b);
322  print (" Y = %8.3e\n",Y);
323 }
Y
#define Y
Definition: radiat.modified_by_Bhargav.c:4
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
menu.c
tuple c
Definition: menu.py:375

Here is the call graph for this function:

int QuarticSolve ( double  b,
double  c,
double  d,
double  e,
double *  z 
)

Solve a quartic equation in the form

\[ z^4 + bz^3 + cz^2 + dz + e = 0 \]

For its purpose, it is assumed that ALL roots are double. This makes things faster.

Parameters
[in]bcoefficient of the quartic
[in]ccoefficient of the quartic
[in]dcoefficient of the quartic
[in]ecoefficient of the quartic
[out]zvector containing the (double) roots of the quartic

Reference:

http://www.1728.com/quartic2.htm

Returns
Return 0 on success, 1 if cubic solver fails, 2 if NaN has been found and 3 if quartic is not satisfied.

Definition at line 23 of file quartic.c.

47 {
48  int status;
49  status = QuarticSolveNew(b,c,d,e,z);
50 
51  if (status != 0){
52  debug_print = 1;
53  QuarticSolveNew(b,c,d,e,z);
54  QUIT_PLUTO(1);
55  }
56 }
QuarticSolveNew
int QuarticSolveNew(double b, double c, double d, double e, double *z)
Definition: quartic.c:175
menu.c
tuple c
Definition: menu.py:375
debug_print
static int debug_print
Definition: quartic.c:20
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call graph for this function:

int QuarticSolveNew ( double  b,
double  c,
double  d,
double  e,
double *  z 
)

Solve a quartic equation in the form

\[ z^4 + bz^3 + cz^2 + dz + e = 0 \]

For its purpose, it is assumed that ALL roots are double. This makes things faster.

Parameters
[in]bcoefficient of the quartic
[in]ccoefficient of the quartic
[in]dcoefficient of the quartic
[in]ecoefficient of the quartic
[out]zvector containing the (double) roots of the quartic

Reference:

https://en.wikipedia.org/wiki/Quartic_function

Returns
Return 0 on success.

Definition at line 175 of file quartic.c.

180 {
181  int n, j, ifail;
182  double b2, sq;
183  double a2, a1, a0, u[4];
184  double p, q, r, f;
185  const double three_256 = 3.0/256.0;
186  const double one_64 = 1.0/64.0;
187  double sz1, sz2, sz3, sz4;
188 
189 
190  b2 = b*b;
191 
192 /* --------------------------------------------------------------
193  1) Compute cubic coefficients using the method outlined in
194  http://eqworld.ipmnet.ru/En/solutions/ae/ae0108.pdf
195  -------------------------------------------------------------- */
196 
197  p = c - b2*0.375;
198  q = d + b2*b*0.125 - b*c*0.5;
199  r = e - 3.0*b2*b2/256.0 + b2*c/16.0 - 0.25*b*d;
200 
201  a2 = 2.0*p;
202  a1 = p*p - 4.0*r;
203  a0 = -q*q;
204 
205  ifail = CubicSolve(a2, a1, a0, u);
206  if (ifail != 0) return 1;
207 
208  u[0] = MAX(u[0],0.0);
209  u[1] = MAX(u[1],0.0);
210  u[2] = MAX(u[2],0.0);
211 
212 
213  if (u[0] != u[0] || u[1] != u[1] || u[2] != u[2]) return 1;
214 
215  sq = -0.5*DSIGN(q);
216  sz1 = sq*sqrt(u[0]);
217  sz2 = sq*sqrt(u[1]);
218  sz3 = sq*sqrt(u[2]);
219 
220  z[0] = -0.25*b + sz1 + sz2 + sz3;
221  z[1] = -0.25*b + sz1 - sz2 - sz3;
222  z[2] = -0.25*b - sz1 + sz2 - sz3;
223  z[3] = -0.25*b - sz1 - sz2 + sz3;
224  SortArray(z,4);
225 
226  if (debug_print){
227  print ("Quartic roots = %f %f %f %f; q = %8.3e\n",z[0],z[1],z[2],z[3],q);
228  CheckSolutions(b,c,d,e,z);
229  }
230 
231  return 0;
232 }
MAX
#define MAX(a, b)
Definition: macros.h:101
CheckSolutions
void CheckSolutions(double b, double c, double d, double e, double *z)
Definition: quartic.c:235
n
static int n
Definition: analysis.c:3
SortArray
void SortArray(double *z, int n)
Definition: quartic.c:326
DSIGN
#define DSIGN(x)
Definition: macros.h:110
CubicSolve
int CubicSolve(double b, double c, double d, double z[])
Definition: quartic.c:59
j
int j
Definition: analysis.c:2
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
menu.c
tuple c
Definition: menu.py:375
debug_print
static int debug_print
Definition: quartic.c:20
q
static Runtime q
Definition: runtime_setup.c:490

Here is the call graph for this function:

Here is the caller graph for this function:

void SortArray ( double *  z,
int  n 
)

Definition at line 326 of file quartic.c.

330 {
331  int j;
332  double f;
333 
334  for (j = 1; j < 4; j++){
335  f = z[j];
336  n = j - 1;
337  while(n >= 0 && z[n] > f) {
338  z[n+1] = z[n];
339  n--;
340  }
341  z[n+1] = f;
342  }
343 }
n
static int n
Definition: analysis.c:3
j
int j
Definition: analysis.c:2

Here is the caller graph for this function:

Variable Documentation

int debug_print =0
static

Definition at line 20 of file quartic.c.