Ice Sheet System Model  4.18
Code documentation
Functions
numerics.h File Reference

prototypes for numerics.h More...

#include "./Verbosity.h"
#include "./GaussPoints.h"
#include "./isnan.h"
#include "./recast.h"
#include "./types.h"
#include "./constants.h"
#include "./OptPars.h"

Go to the source code of this file.

Functions

IssmDouble min (IssmDouble a, IssmDouble b)
 
IssmDouble max (IssmDouble a, IssmDouble b)
 
int min (int a, int b)
 
int max (int a, int b)
 
void BrentSearch (IssmDouble **pJ, OptPars optpars, IssmDouble *X0, IssmDouble(*f)(IssmDouble *, void *), IssmDouble(*g)(IssmDouble **, IssmDouble *, void *), void *usr)
 
void cross (IssmDouble *result, IssmDouble *vector1, IssmDouble *vector2)
 
void XZvectorsToCoordinateSystem (IssmDouble *T, IssmDouble *xzvectors)
 
int cubic (IssmDouble a, IssmDouble b, IssmDouble c, IssmDouble d, IssmDouble X[3], int *num)
 
IssmDouble legendre (IssmDouble Pn1, IssmDouble Pn2, IssmDouble x, int n)
 
IssmDouble * p_polynomial_value (int m, int n, IssmDouble *x)
 
int NewtonSolveDnorm (IssmDouble *pdnorm, IssmDouble c1, IssmDouble c2, IssmDouble c3, IssmDouble n, IssmDouble dnorm)
 
IssmDouble ODE1 (IssmDouble alpha, IssmDouble beta, IssmDouble Si, IssmDouble dt, int method)
 

Detailed Description

prototypes for numerics.h

Definition in file numerics.h.

Function Documentation

◆ min() [1/2]

IssmDouble min ( IssmDouble  a,
IssmDouble  b 
)

Definition at line 14 of file extrema.cpp.

14  {
15  if (a<b)return a;
16  else return b;
17 }

◆ max() [1/2]

IssmDouble max ( IssmDouble  a,
IssmDouble  b 
)

Definition at line 24 of file extrema.cpp.

24  {
25  if (a>b)return a;
26  else return b;
27 }

◆ min() [2/2]

int min ( int  a,
int  b 
)

Definition at line 19 of file extrema.cpp.

19  {
20  if (a<b)return a;
21  else return b;
22 }

◆ max() [2/2]

int max ( int  a,
int  b 
)

Definition at line 29 of file extrema.cpp.

29  {
30  if (a>b)return a;
31  else return b;
32 }

◆ BrentSearch()

void BrentSearch ( IssmDouble **  pJ,
OptPars  optpars,
IssmDouble *  X0,
IssmDouble(*)(IssmDouble *, void *)  f,
IssmDouble(*)(IssmDouble **, IssmDouble *, void *)  g,
void *  usr 
)

Definition at line 23 of file BrentSearch.cpp.

23  {
24 
25  /* This routine is optimizing a given function using Brent's method
26  * (Golden or parabolic procedure)*/
27 
28  /*Intermediary*/
29  int iter;
30  IssmDouble si,gold,intervalgold,oldintervalgold;
31  IssmDouble parab_num,parab_den,distance;
32  IssmDouble fxmax,fxmin,fxbest;
33  IssmDouble fx,fx1,fx2;
34  IssmDouble x,x1,x2,xm,xbest;
35  IssmDouble tol1,tol2,seps;
36  IssmDouble tolerance = 1.e-4;
37 
38  /*Recover parameters:*/
39  int nsteps = optpars.nsteps;
40  int nsize = optpars.nsize;
41  IssmDouble xmin = optpars.xmin;
42  IssmDouble xmax = optpars.xmax;
43  int *maxiter = optpars.maxiter;
44  IssmDouble *cm_jump = optpars.cm_jump;
45 
46  /*Initialize gradient and controls*/
47  IssmDouble* G = NULL;
48  IssmDouble* J = xNew<IssmDouble>(nsteps);
49  IssmDouble* X = xNew<IssmDouble>(nsize);
50 
51  /*Header of printf*/
52  _printf0_("\n");
53  _printf0_(" x | Cost function f(x) | List of contributions\n");
54 
55  /*start iterations*/
56  for(int n=0;n<nsteps;n++){
57 
58  /*Print iteration number*/
59  _printf0_("====================== step "<< n+1 << "/" << nsteps <<" ===============================\n");
60 
61  /*Reset some variables*/
62  iter = 0;
63  xmin = 0.;
64  xmax = 1.;
65  bool loop = true;
66  cout<<setprecision(5);
67 
68  /*Get current Gradient at xmin=0*/
69  _printf0_(" x = "<<setw(9)<<xmin<<" | ");
70  fxmin = (*g)(&G,X0,usr); if(xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
71 
72  /*Get f(xmax)*/
73  _printf0_(" x = "<<setw(9)<<xmax<<" | ");
74  for(int i=0;i<nsize;i++) X[i]=X0[i]+xmax*G[i];
75  fxmax = (*f)(X,usr); if (xIsNan<IssmDouble>(fxmax)) _error_("Function evaluation returned NaN");
76  //if(VerboseControl()) _printf0_(" N/A "<<setw(12)<<xmax<<" "<<setw(12)<<fxmax<<" N/A boundary\n");
77 
78  /*test if jump option activated and xmin==0*/
79  if(!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && (fxmax/fxmin)<cm_jump[n]){
80  for(int i=0;i<nsize;i++) X0[i]=X0[i]+xmax*G[i];
81  xDelete<IssmDouble>(G);
82  J[n]=fxmax;
83  continue;
84  }
85 
86  /*initialize optimization variables*/
87  seps=sqrt(DBL_EPSILON); //precision of a IssmDouble
88  distance=0.0; //new_x=old_x + distance
89  gold=0.5*(3.0-sqrt(5.0)); //gold = 1 - golden ratio
90  intervalgold=0.0; //distance used by Golden procedure
91 
92  /*1: initialize the values of the 4 x needed (x1,x2,x,xbest)*/
93  x1=xmin+gold*(xmax-xmin);
94  x2=x1;
95  xbest=x1;
96  x=xbest;
97 
98  /*2: call the function to be evaluated*/
99  _printf0_(" x = "<<setw(9)<<x<<" | ");
100  for(int i=0;i<nsize;i++) X[i]=X0[i]+x*G[i];
101  fxbest = (*f)(X,usr); if(xIsNan<IssmDouble>(fxbest)) _error_("Function evaluation returned NaN");
102  iter++;
103 
104  /*3: update the other variables*/
105  fx1=fxbest;
106  fx2=fxbest;
107  /*xm is always in the middle of a and b*/
108  xm=0.5*(xmin+xmax);
109  /*update tolerances*/
110  tol1=seps*sqrt(pow(xbest,2))+tolerance/3.0;
111  tol2=2.0*tol1;
112 
113  /*4: print result*/
114  //if(VerboseControl())
115  //_printf0_(" "<<setw(5)<<iter<<" "<<setw(12)<<xbest<<" "<<setw(12)<<fxbest<<" "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<" initial\n");
116  if (!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && ((fxbest/fxmin)<cm_jump[n])){
117  //if(VerboseControl()) _printf0_(" optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump[n] << "\n");
118  loop=false;
119  }
120 
121  while(loop){
122 
123  bool goldenflag=true;
124 
125  // Is a parabolic fit possible ?
126  if (sqrt(pow(intervalgold,2))>tol1){
127 
128  // Yes, so fit parabola
129  goldenflag=false;
130  parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);;
131  parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1);
132 
133  //reverse p if necessary
134  if(parab_den>0.0){
135  parab_num=-parab_num;
136  }
137  parab_den=sqrt(pow(parab_den,2));
138  oldintervalgold=intervalgold;
139  intervalgold=distance;
140 
141  // Is the parabola acceptable (we use seps here because in some configuration parab_num==parab_den*(xmax-xbest)
142  // and the result is not repeatable anymore
143  if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) &&
144  (parab_num>parab_den*(xmin-xbest)+seps) &&
145  (parab_num<parab_den*(xmax-xbest)-seps)){
146 
147  // Yes, parabolic interpolation step
148  distance=parab_num/parab_den;
149  x=xbest+distance;
150 
151  // f must not be evaluated too close to min_x or max_x
152  if (((x-xmin)<tol2) || ((xmax-x)<tol2)){
153  if ((xm-xbest)<0.0) si=-1;
154  else si=1;
155  //compute new distance
156  distance=tol1*si;
157  }
158  }
159  else{
160  // Not acceptable, must do a golden section step
161  goldenflag=true;
162  }
163  }
164 
165  //Golden procedure
166  if(goldenflag){
167  // compute the new distance d
168  if(xbest>=xm){
169  intervalgold=xmin-xbest;
170  }
171  else{
172  intervalgold=xmax-xbest;
173  }
174  distance=gold*intervalgold;
175  }
176 
177  // The function must not be evaluated too close to xbest
178  if(distance<0) si=-1;
179  else si=1;
180  if(sqrt(pow(distance,2))>tol1) x=xbest+si*sqrt(pow(distance,2));
181  else x=xbest+si*tol1;
182 
183  //evaluate function on x
184  _printf0_(" x = "<<setw(9)<<x<<" | ");
185  for(int i=0;i<nsize;i++) X[i]=X0[i]+x*G[i];
186  fx = (*f)(X,usr); if(xIsNan<IssmDouble>(fx)) _error_("Function evaluation returned NaN");
187  iter++;
188 
189  // Update a, b, xm, x1, x2, tol1, tol2
190  if (fx<=fxbest){
191  if (x>=xbest) xmin=xbest;
192  else xmax=xbest;
193  x1=x2; fx1=fx2;
194  x2=xbest; fx2=fxbest;
195  xbest=x; fxbest=fx;
196  }
197  else{ // fx > fxbest
198  if (x<xbest) xmin=x;
199  else xmax=x;
200  if ((fx<=fx2) || (x2==xbest)){
201  x1=x2; fx1=fx2;
202  x2=x; fx2=fx;
203  }
204  else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){
205  x1=x; fx1=fx;
206  }
207  }
208  xm = 0.5*(xmin+xmax);
209  tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0;
210  tol2=2.0*tol1;
211  //if(VerboseControl())
212  // _printf0_(" "<<setw(5)<<iter<<" "<<setw(12)<<x<<" "<<setw(12)<<fx<<" "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<
213  // " "<<(goldenflag?"golden":"parabolic")<<"\n");
214 
215  /*Stop the optimization?*/
216  if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(xmax-xmin))){
217  //if(VerboseControl()) _printf0_(" optimization terminated: current x satisfies criteria 'tolx'=" << tolerance << "\n");
218  loop=false;
219  }
220  else if (iter>=maxiter[n]){
221  //if(VerboseControl()) _printf0_(" exiting: Maximum number of iterations has been exceeded ('maxiter'=" << maxiter[n] << ")\n");
222  loop=false;
223  }
224  else if (!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && ((fxbest/fxmin)<cm_jump[n])){
225  //if(VerboseControl()) _printf0_(" optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump[n] << "\n");
226  loop=false;
227  }
228  else{
229  //continue
230  loop=true;
231  }
232  }//end while
233 
234  //Now, check that the value on the boundaries are not better than current fxbest
235  if (fxbest>fxmin){
236  xbest=optpars.xmin; fxbest=fxmin;
237  }
238  if (fxbest>fxmax){
239  xbest=optpars.xmax; fxbest=fxmax;
240  }
241 
242  /*Assign output pointers: */
243  for(int i=0;i<nsize;i++) X0[i]=X0[i]+xbest*G[i];
244  xDelete<IssmDouble>(G);
245  J[n]=fxbest;
246  }
247 
248  /*return*/
249  xDelete<IssmDouble>(X);
250  *pJ=J;
251 }

◆ cross()

void cross ( IssmDouble *  result,
IssmDouble *  vector1,
IssmDouble *  vector2 
)

Definition at line 13 of file cross.cpp.

13  {
14 
15  /*result,vector1 and vector2 are all assumed to be of size 3: */
16 
17  result[0]=vector1[1]*vector2[2]-vector1[2]*vector2[1];
18  result[1]=vector1[2]*vector2[0]-vector1[0]*vector2[2];
19  result[2]=vector1[0]*vector2[1]-vector1[1]*vector2[0];
20 
21 }

◆ XZvectorsToCoordinateSystem()

void XZvectorsToCoordinateSystem ( IssmDouble *  T,
IssmDouble *  xzvectors 
)

Definition at line 8 of file XZvectorsToCoordinateSystem.cpp.

8  {
9 
10  IssmDouble x[3],y[3],z[3];
11  IssmDouble x_norm, y_norm, z_norm;
12 
13  for(int i=0;i<6;i++){
14  if(xIsNan<IssmDouble>(xzvectors[i])){
15  /*At least one NaN found: default to Id*/
16  T[0*3+0] = 1.0; T[0*3+1] = 0.0; T[0*3+2] = 0.0;
17  T[1*3+0] = 0.0; T[1*3+1] = 1.0; T[1*3+2] = 0.0;
18  T[2*3+0] = 0.0; T[2*3+1] = 0.0; T[2*3+2] = 1.0;
19 
20  return;
21  }
22  }
23 
24  /* get input {x} (vector in local x-z plane): */
25  x[0] = xzvectors[0];
26  x[1] = xzvectors[1];
27  x[2] = xzvectors[2];
28 
29  /* get input {z} (local tangent plane normal vector): */
30  z[0] = xzvectors[3];
31  z[1] = xzvectors[4];
32  z[2] = xzvectors[5];
33 
34  /* compute {y} = {z} x {x}: */
35  y[0] = x[2]*z[1] - x[1]*z[2];
36  y[1] = -x[2]*z[0] + x[0]*z[2];
37  y[2] = x[1]*z[0] - x[0]*z[1];
38 
39  /* normalise {x}, {y} and {z} to form unit vectors {i_hat}, {j_hat} and {k_hat};
40  store in {x}, {y}, and {z}: */
41  x_norm = sqrt( x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
42  y_norm = sqrt( y[0]*y[0] + y[1]*y[1] + y[2]*y[2]);
43  z_norm = sqrt( z[0]*z[0] + z[1]*z[1] + z[2]*z[2]);
44 
45  x[0] = x[0]/x_norm; x[1] = x[1]/x_norm; x[2] = x[2]/x_norm;
46  y[0] = y[0]/y_norm; y[1] = y[1]/y_norm; y[2] = y[2]/y_norm;
47  z[0] = z[0]/z_norm; z[1] = z[1]/z_norm; z[2] = z[2]/z_norm;
48 
49  /* Tlg columns are just {i_hat}, {j_hat} and {k_hat}, respectively: */
50  T[0*3+0] = x[0]; T[0*3+1] = y[0]; T[0*3+2] = z[0];
51  T[1*3+0] = x[1]; T[1*3+1] = y[1]; T[1*3+2] = z[1];
52  T[2*3+0] = x[2]; T[2*3+1] = y[2]; T[2*3+2] = z[2];
53 }

◆ cubic()

int cubic ( IssmDouble  a,
IssmDouble  b,
IssmDouble  c,
IssmDouble  d,
IssmDouble  X[3],
int *  num 
)

Definition at line 21 of file cubic.cpp.

21  {
22  /* Find the real roots of linear/quadratic and cubic equations:
23  *
24  * a x^3 + bx^2 + cx + d = 0
25  *
26  * returns the roots in x
27  * num is the number of roots */
28 
29  /*Some useful constants*/
30  const IssmDouble pi = 3.1415926535897932;
31  const IssmDouble third = 1./3.;
32 
33  /*Intermediaries*/
34  IssmDouble U[3],W,P,Q,delta,phi;
35 
36  /* determine the degree of the polynomial */
37  if (a != 0.0){
38  //cubic problem
39  W = b/a *third;
40  P = pow((c/a *third - pow(W,2)),3);
41  Q = -.5 *(2.0*pow(W,3)-(c*W-d)/a );
42  delta = pow(Q,2)+P;
43  if ( delta < 0.0 ){
44  //three real solutions!
45  //Confine the argument of coeffCOS to the interval [-1;1]!
46  phi = acos(min(1.0,max(-1.0,Q/sqrt(-P))));
47  P = 2.0*pow((-P),(5.e-1*third));
48  for(int i=0;i<3;i++) U[i] = P*cos((phi+2*((IssmDouble)i)*pi)*third)-W;
49  x[0] = min(U[0], min(U[1], U[2]));
50  x[1] = max(min(U[0], U[1]),max( min(U[0], U[2]), min(U[1], U[2])));
51  x[2] = max(U[0], max(U[1], U[2]));
52  *num = 3;
53  }
54  else{
55  // only one real solution!
56  delta = sqrt(delta);
57  x[0] = CBRT(Q+delta)+CBRT(Q-delta)-W;
58  *num=1;
59  }
60  }
61  else if (b != 0.0){
62  // quadratic problem
63  P = 0.5*c/b;
64  delta = pow(P,2)-d/b;
65  if (delta > 0.0){
66  // 2 real solutions
67  x[0] = -P - sqrt(delta);
68  x[1] = -P + sqrt(delta);
69  *num = 2;
70  }
71  else{
72  // no real solution
73  *num = 0;
74  }
75  }
76  else if (c != 0.0){
77  //linear equation
78  x[0] = d/c;
79  *num = 1;
80  }
81  else{
82  //no equation
83  *num = 0;
84  }
85 
86  /* perform one step of a newton iteration in order to minimize round-off errors */
87  for(int i=0;i<*num;i++){
88  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));
89  }
90 
91  return 0;
92 }

◆ legendre()

IssmDouble legendre ( IssmDouble  Pn1,
IssmDouble  Pn2,
IssmDouble  x,
int  n 
)

◆ p_polynomial_value()

IssmDouble* p_polynomial_value ( int  m,
int  n,
IssmDouble *  x 
)

Definition at line 142 of file legendre.cpp.

142  {
143 
144  /******************************************************************************{{{/
145  Purpose:
146 
147  P_POLYNOMIAL_VALUE evaluates the Legendre polynomials P(n,x).
148 
149  Discussion:
150 
151  P(n,1) = 1.
152  P(n,-1) = (-1)^N.
153  | P(n,x) | <= 1 in [-1,1].
154 
155  The N zeroes of P(n,x) are the abscissas used for Gauss-Legendre
156  quadrature of the integral of a function F(X) with weight function 1
157  over the interval [-1,1].
158 
159  The Legendre polynomials are orthogonal under the inner product defined
160  as integration from -1 to 1:
161 
162  Integral ( -1 <= X <= 1 ) P(I,X) * P(J,X) dX
163  = 0 if I =/= J
164  = 2 / ( 2*I+1 ) if I = J.
165 
166  Except for P(0,X), the integral of P(I,X) from -1 to 1 is 0.
167 
168  A function F(X) defined on [-1,1] may be approximated by the series
169  C0*P(0,x) + C1*P(1,x) + ... + CN*P(n,x)
170  where
171  C(I) = (2*I+1)/(2) * Integral ( -1 <= X <= 1 ) F(X) P(I,x) dx.
172 
173  The formula is:
174 
175  P(n,x) = (1/2^N) * sum ( 0 <= M <= N/2 ) C(N,M) C(2N-2M,N) X^(N-2*M)
176 
177  Differential equation:
178 
179  (1-X*X) * P(n,x)'' - 2 * X * P(n,x)' + N * (N+1) = 0
180 
181  First terms:
182 
183  P( 0,x) = 1
184  P( 1,x) = 1 X
185  P( 2,x) = ( 3 X^2 - 1)/2
186  P( 3,x) = ( 5 X^3 - 3 X)/2
187  P( 4,x) = ( 35 X^4 - 30 X^2 + 3)/8
188  P( 5,x) = ( 63 X^5 - 70 X^3 + 15 X)/8
189  P( 6,x) = ( 231 X^6 - 315 X^4 + 105 X^2 - 5)/16
190  P( 7,x) = ( 429 X^7 - 693 X^5 + 315 X^3 - 35 X)/16
191  P( 8,x) = ( 6435 X^8 - 12012 X^6 + 6930 X^4 - 1260 X^2 + 35)/128
192  P( 9,x) = (12155 X^9 - 25740 X^7 + 18018 X^5 - 4620 X^3 + 315 X)/128
193  P(10,x) = (46189 X^10-109395 X^8 + 90090 X^6 - 30030 X^4 + 3465 X^2-63)/256
194 
195  Recursion:
196 
197  P(0,x) = 1
198  P(1,x) = x
199  P(n,x) = ( (2*n-1)*x*P(n-1,x)-(n-1)*P(n-2,x) ) / n
200 
201  P'(0,x) = 0
202  P'(1,x) = 1
203  P'(N,x) = ( (2*N-1)*(P(N-1,x)+X*P'(N-1,x)-(N-1)*P'(N-2,x) ) / N
204 
205  Licensing:
206 
207  This code is distributed under the GNU LGPL license.
208 
209  Modified:
210 
211  08 August 2013
212 
213  Author:
214 
215  John Burkardt
216 
217  Reference:
218 
219  Milton Abramowitz, Irene Stegun,
220  Handbook of Mathematical Functions,
221  National Bureau of Standards, 1964,
222  ISBN: 0-486-61272-4,
223  LC: QA47.A34.
224 
225  Daniel Zwillinger, editor,
226  CRC Standard Mathematical Tables and Formulae,
227  30th Edition,
228  CRC Press, 1996.
229 
230  Parameters:
231 
232  Input, int M, the number of evaluation points.
233 
234  Input, int N, the highest order polynomial to evaluate.
235  Note that polynomials 0 through N will be evaluated.
236 
237  Input, double X[M], the evaluation points.
238 
239  Output, double P_POLYNOMIAL_VALUE[M*(N+1)], the values of the Legendre
240  polynomials of order 0 through N.
241  }}}*/
242 
243  int i,j;
244 
245  if(n<0) return NULL;
246 
247  IssmDouble* v = xNew<IssmDouble>(m*(n+1));
248 
249  for ( i = 0; i < m; i++ ) v[i*(n+1)+0] = 1.0;
250  if ( n < 1 ) return v;
251 
252  for ( i = 0; i < m; i++ ) v[i*(n+1)+1] = x[i];
253 
254  for ( i = 0; i < m; i++ ) {
255  for ( j = 2; j <= n; j++ ) {
256 
257  v[j+i*(n+1)] = ( ( IssmDouble ) ( 2 * j - 1 ) * x[i] * v[(j-1)+i*(n+1)]
258  - ( IssmDouble ) ( j - 1 ) * v[(j-2)+i*(n+1)] )
259  / ( IssmDouble ) ( j );
260 
261  }
262  }
263 
264  return v;
265 }

◆ NewtonSolveDnorm()

int NewtonSolveDnorm ( IssmDouble *  pdnorm,
IssmDouble  c1,
IssmDouble  c2,
IssmDouble  c3,
IssmDouble  n,
IssmDouble  dnorm 
)

Definition at line 5 of file NewtonSolveDnorm.cpp.

5  {
6  /* solve the following equation using Newton-Raphson
7  *
8  * c1*x^(s-1) + c2*x = c3
9  *
10  * s = (1+n)/n
11  *
12  * we solve y = 10^x:
13  */
14 
15  /*trivial solution*/
16  if(c3==0.){
17  *pdnorm = 0.;
18  return 0;
19  }
20 
21  /*Intermediaries*/
22  int counter = 0;
23  IssmDouble s = (1.+n)/n;
24  IssmDouble y2;
25  IssmDouble threshold = 1.e-12;
26 
27  /*Initial guess*/
28  _assert_(dnorm>0.);
29  IssmDouble y1 = log10(dnorm);
30 
31  while(true){
32 
33  /*Newton step*/
34  y2 = y1 - (c1*pow(pow(10.,y1),s-1.) + c2*pow(10.,y1) - c3)/((s-1)*c1*log(10.)*pow(pow(10.,y1),s-1.) + c2*log(10.)*pow(10.,y1));
35 
36  if( fabs(y2-y1)/(fabs(y2))<threshold ){
37  break;
38  }
39  else{
40  y1=y2;
41  counter++;
42  }
43 
44  if(counter>50) break;
45  }
46 
47  /*Avoid extremely large values that indicate non convergence*/
48  if(y2>50.) y2 = 50;
49 
50  /*Assign output pointer*/
51  *pdnorm = pow(10.,y2);
52  return 0;
53 }

◆ ODE1()

IssmDouble ODE1 ( IssmDouble  alpha,
IssmDouble  beta,
IssmDouble  Si,
IssmDouble  dt,
int  method 
)

Definition at line 5 of file ODE1.cpp.

5  {
6  /* solve the following equation:
7  *
8  * dS/dt = alpha S + beta
9  *
10  * method 0: Forward Euler (explicit)
11  * method 1: backward Euler (implicit)
12  * method 2: Crank Nicolson
13  *
14  * return S^{i+1} based on Si, dt, alpha and beta
15  */
16 
17  switch(method){
18  case 0: return Si*(1.+alpha*dt) + beta*dt;
19  case 1: return (Si+beta*dt)/(1.-alpha*dt);
20  case 2: return (Si*(1.+alpha*dt/2.) + beta*dt)/(1-alpha*dt/2.);
21  default: _error_("not supported yet");
22  }
23 }
OptPars::nsize
int nsize
Definition: OptPars.h:17
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
OptPars::maxiter
int * maxiter
Definition: OptPars.h:15
OptPars::xmax
IssmDouble xmax
Definition: OptPars.h:13
fx
IssmDouble fx(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
Definition: fsanalyticals.cpp:218
alpha
IssmDouble alpha(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
Definition: fsanalyticals.cpp:221
OptPars::nsteps
int nsteps
Definition: OptPars.h:16
OptPars::xmin
IssmDouble xmin
Definition: OptPars.h:12
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
min
IssmDouble min(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:14
OptPars::cm_jump
IssmDouble * cm_jump
Definition: OptPars.h:14
CBRT
IssmDouble CBRT(IssmDouble Z)
Definition: cubic.cpp:5
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
Generated on Thu Jul 2 2020 08:09:25 for Ice Sheet System Model by   doxygen 1.8.19