Changeset 13061
- Timestamp:
- 08/16/12 09:24:06 (13 years ago)
- Location:
- issm/trunk-jpl/src/c/shared/Numerics
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/shared/Numerics/cubic.cpp
r13050 r13061 2 2 #include "./numerics.h" 3 3 4 int signR(double Z){ 5 int ret; 6 if (Z > 0.0) ret = 1; 7 if (Z < 0.0) ret = -1; 8 if (Z == 0.0) ret =0; 4 int cubic(IssmDouble a,IssmDouble b,IssmDouble c,IssmDouble d, IssmDouble x[3], int* num){ 5 /* Find the real roots of linear/quadratic and cubic equations: 6 * 7 * a x^3 + bx^2 + cx + d = 0 8 * 9 * returns the roots in x 10 * num is the number of roots */ 9 11 10 return ret; 11 } 12 /*Some useful constants*/ 13 const IssmDouble pi = 3.1415926535897932; 14 const IssmDouble third = 1./3.; 12 15 13 double CBRT(double Z){ 14 double ret; 15 const double THIRD = 1./3.; 16 //define cubic root as statement function 17 //SIGN has different meanings in both C and Fortran 18 // Was unable to use the sign command of C, so wrote my own 19 // that why a new variable needs to be introduced that keeps track of the sign of 20 // SIGN is supposed to return a 1, -1 or 0 depending on what the sign of the argument is 21 ret = fabs(pow(fabs(Z),THIRD)) * signR(Z); 22 return ret; 23 } 16 /*Intermediaries*/ 17 IssmDouble U[3],W,P,Q,delta,phi; 24 18 25 int cubic(double A[4], double X[3], int* num){ 26 const double PI = 3.1415926535897932; 27 const double THIRD = 1./3.; 28 double U[3],W, P, Q, DIS, PHI; 29 int i; 30 31 //define cubic root as statement function 32 // In C, the function is defined outside of the cubic fct 33 34 // ====determine the degree of the polynomial ==== 35 36 if (A[3] != 0.0){ 19 /* determine the degree of the polynomial */ 20 if (a != 0.0){ 37 21 //cubic problem 38 W = A[2]/A[3]*THIRD;39 P = pow((A[1]/A[3]*THIRD- pow(W,2)),3);40 Q = -.5*(2.0*pow(W,3)-(A[1]*W-A[0])/A[3]);41 DIS= pow(Q,2)+P;42 if ( DIS< 0.0 ){22 W = b/a *third; 23 P = pow((c/a *third - pow(W,2)),3); 24 Q = -.5 *(2.0*pow(W,3)-(c*W-d)/a ); 25 delta = pow(Q,2)+P; 26 if ( delta < 0.0 ){ 43 27 //three real solutions! 44 //Confine the argument of ACOS to the interval [-1;1]!45 PHI= acos(min(1.0,max(-1.0,Q/sqrt(-P))));46 P =2.0*pow((-P),(5.e-1*THIRD));47 for (i=0;i<3;i++) U[i] = P*cos((PHI+2*((double)i)*PI)*THIRD)-W;48 X[0] = min(U[0], min(U[1], U[2]));49 X[1] = max(min(U[0], U[1]),max( min(U[0], U[2]), min(U[1], U[2])));50 X[2] = max(U[0], max(U[1], U[2]));28 //Confine the argument of coeffCOS to the interval [-1;1]! 29 phi = acos(min(1.0,max(-1.0,Q/sqrt(-P)))); 30 P = 2.0*pow((-P),(5.e-1*third)); 31 for(int i=0;i<3;i++) U[i] = P*cos((phi+2*((IssmDouble)i)*pi)*third)-W; 32 x[0] = min(U[0], min(U[1], U[2])); 33 x[1] = max(min(U[0], U[1]),max( min(U[0], U[2]), min(U[1], U[2]))); 34 x[2] = max(U[0], max(U[1], U[2])); 51 35 *num = 3; 52 36 } 53 37 else{ 54 38 // only one real solution! 55 DIS = sqrt(DIS);56 X[0] = CBRT(Q+DIS)+CBRT(Q-DIS)-W;39 delta = sqrt(delta); 40 x[0] = pow(Q+delta,third)+pow(Q-delta,third)-W; 57 41 *num=1; 58 42 } 59 43 } 60 else if ( A[2]!= 0.0){44 else if (b != 0.0){ 61 45 // quadratic problem 62 P = 0.5*A[1]/A[2]; 63 DIS = pow(P,2)-A[0]/A[2]; 64 if (DIS > 0.0) 65 { 46 P = 0.5*c/b; 47 delta = pow(P,2)-d/b; 48 if (delta > 0.0){ 66 49 // 2 real solutions 67 X[0] = -P - sqrt(DIS);68 X[1] = -P + sqrt(DIS);69 *num =2;50 x[0] = -P - sqrt(delta); 51 x[1] = -P + sqrt(delta); 52 *num = 2; 70 53 } 71 else 72 { 54 else{ 73 55 // no real solution 74 *num =0;56 *num = 0; 75 57 } 76 58 } 77 else if (A[1] != 0.0) 78 { 59 else if (c != 0.0){ 79 60 //linear equation 80 X[0] =A[0]/A[1];81 *num =1;61 x[0] = d/c; 62 *num = 1; 82 63 } 83 else 84 { 64 else{ 85 65 //no equation 86 *num =0;66 *num = 0; 87 67 } 88 /* 89 * ==== perform one step of a newton iteration in order to minimize 90 * round-off errors ==== 91 */ 92 for (i=0;i<*num;i++){ 93 X[i] = X[i] - (A[0]+X[i]*(A[1]+X[i]*(A[2]+X[i]*A[3])))/(A[1]+X[i]*(2.0*A[2]+X[i]*3.0*A[3])); 68 69 /* perform one step of a newton iteration in order to minimize round-off errors */ 70 for(int i=0;i<*num;i++){ 71 x[i] = x[i] - (d+x[i]*(c+x[i]*(b+x[i]*a)))/(c+x[i]*(2.0*b+x[i]*3.0*a)); 94 72 } 73 95 74 return 0; 96 75 } -
issm/trunk-jpl/src/c/shared/Numerics/numerics.h
r13050 r13061 17 17 struct OptPars; 18 18 19 IssmDouble min(IssmDouble a,IssmDouble b);20 IssmDouble max(IssmDouble a,IssmDouble b);21 int min(int a,int b);22 int max(int a,int b);23 IssmDouble OptFunc(IssmDouble scalar, OptArgs*optargs);24 void BrentSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs*optargs);25 void OptimalSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs*optargs);26 void cross(IssmDouble* result,IssmDouble* vector1,IssmDouble*vector2);27 void IsInputConverged(IssmDouble* peps, Input** new_inputs,Input**old_inputs,int num_inputs,int criterion_enum);28 void UnitConversion(IssmDouble*values, int numvalues,int direction_enum, int type_enum);29 IssmDouble UnitConversion(IssmDouble value, int direction_enum, int type_enum);30 char * OptionsFromAnalysis(Parameters*parameters,int analysis_type);31 void XZvectorsToCoordinateSystem(IssmDouble* T,IssmDouble*xzvectors);32 int cubic(double A[4], double X[3], int*num);19 IssmDouble min(IssmDouble a,IssmDouble b); 20 IssmDouble max(IssmDouble a,IssmDouble b); 21 int min(int a,int b); 22 int max(int a,int b); 23 IssmDouble OptFunc(IssmDouble scalar, OptArgs *optargs); 24 void BrentSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs*optargs); 25 void OptimalSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs*optargs); 26 void cross(IssmDouble *result,IssmDouble*vector1,IssmDouble*vector2); 27 void IsInputConverged(IssmDouble *peps, Input**new_inputs,Input**old_inputs,int num_inputs,int criterion_enum); 28 void UnitConversion(IssmDouble *values, int numvalues,int direction_enum, int type_enum); 29 IssmDouble UnitConversion(IssmDouble value, int direction_enum, int type_enum); 30 char *OptionsFromAnalysis(Parameters *parameters,int analysis_type); 31 void XZvectorsToCoordinateSystem(IssmDouble *T,IssmDouble*xzvectors); 32 int cubic(IssmDouble a, IssmDouble b, IssmDouble c, IssmDouble d, double X[3], int *num); 33 33 #ifdef _HAVE_PETSC_ 34 34 void PetscOptionsFromAnalysis(Parameters* parameters,int analysis_type);
Note:
See TracChangeset
for help on using the changeset viewer.