NumericalMath


1   package JSci.maths;
2   
3   import JSci.maths.analysis.RealFunction;
4   import JSci.maths.analysis.RealFunction2D;
5   
6   /**
7   * The numerical math library.
8   * This class cannot be subclassed or instantiated because all methods are static.
9   * @version 0.9
10  * @author Mark Hale
11  */
12  public final class NumericalMath extends AbstractMath {
13          private NumericalMath() {}
14  
15          /**
16          * Calculates the roots of the quadratic equation
17          * ax<sup>2</sup>+bx+c=0.
18          * @return an array containing the two roots.
19          */
20          public static double[] solveQuadratic(final double a,final double b,final double c) {
21                  final double roots[]=new double[2];
22                  double q;
23                  if(b < 0.0)
24                          q = (-b+Math.sqrt(b*b-4.0*a*c))/2.0;
25                  else
26                          q = (-b-Math.sqrt(b*b-4.0*a*c))/2.0;
27                  roots[0] = q/a;
28                  roots[1] = c/q;
29                  return roots;
30          }
31          /**
32          * Finds a root using the bisection method.
33          * @param a lower bound.
34          * @param b upper bound.
35          */
36          public static double bisection(Mapping func, double a, double b, final double tol) {
37                  final int signa = ExtraMath.sign(func.map(a));
38                  final int signb = ExtraMath.sign(func.map(b));
39                  if(signa == signb)
40                          throw new IllegalArgumentException  ("Bounds do not bracket a root.");
41                  double x;
42                  do {
43                          x = (a + b)/2.0;
44                          int signx = ExtraMath.sign(func.map(x));
45                          if(signx == signa) {
46                                  a = x;
47                          } else if(signx == signb) {
48                                  b = x;
49                          } else {
50                                  a = x;
51                                  b = x;
52                          }
53                  } while(Math.abs(a-b) > tol);
54                  return x;
55          }
56          /**
57          * Finds a root using the false position method.
58          * @param a lower bound.
59          * @param b upper bound.
60          */
61          public static double falsePosition(Mapping func, double a, double b, final double tol) {
62                  double fa = func.map(a);
63                  double fb = func.map(b);
64                  final int signa = ExtraMath.sign(fa);
65                  final int signb = ExtraMath.sign(fb);
66                  if(signa == signb)
67                          throw new IllegalArgumentException  ("Bounds do not bracket a root.");
68                  double x;
69                  double delta;
70                  do {
71                          x = a - (b-a)*fa/(fb-fa);
72                          double fx = func.map(x);
73                          int signx = ExtraMath.sign(fx);
74                          if(signx == signa) {
75                                  delta = x-a;
76                                  a = x;
77                                  fa = fx;
78                          } else if(signx == signb) {
79                                  delta = b-x;
80                                  b = x;
81                                  fb = fx;
82                          } else {
83                                  delta = 0.0;
84                          }
85                  } while(Math.abs(delta) > tol);
86                  return x;
87          }
88          /**
89          * Finds a root using the Newton-Raphson method.
90          * @param x initial guess.
91          */
92          public static double newtonRaphson(RealFunction func, double x, final int maxIter, final double tol) throws MaximumIterationsExceededException {
93                  RealFunction deriv = func.differentiate();
94                  double delta;
95                  int n = 0;
96                  do {
97                          delta = -func.map(x)/deriv.map(x);
98                          x += delta;
99                          if(++n > maxIter)
100                                 throw new MaximumIterationsExceededException("No convergence after "+maxIter+" iterations.");
101                 } while(Math.abs(delta) > tol);
102                 return x;
103         }
104         /**
105         * Uses the Euler method to solve an ODE.
106         * @param y an array to be filled with y values, set y[0] to initial condition.
107         * @param func dy/dt as a function of y.
108         * @param dt step size.
109         * @return y.
110         */
111         public static double[] euler(final double y[],final Mapping func,final double dt) {
112                 for(int i=0;i<y.length-1;i++)
113                         y[i+1]=y[i]+dt*func.map(y[i]);
114                 return y;
115         }
116         /**
117         * Uses the Leap-Frog method to solve an ODE.
118         * @param y an array to be filled with y values, set y[0], y[1] to initial conditions.
119         * @param func dy/dt as a function of y.
120         * @param dt step size.
121         * @return y.
122         */
123         public static double[] leapFrog(final double y[],final Mapping func,final double dt) {
124                 final double two_dt = 2.0*dt;
125                 for(int i=1;i<y.length-1;i++)
126                         y[i+1]=y[i-1]+two_dt*func.map(y[i]);
127                 return y;
128         }
129         /**
130         * Uses the 2nd order Runge-Kutta method to solve an ODE.
131         * @param y an array to be filled with y values, set y[0] to initial condition.
132         * @param func dy/dt as a function of y.
133         * @param dt step size.
134         * @return y.
135         */
136         public static double[] rungeKutta2(final double y[],final Mapping func,final double dt) {
137                 final double dt2 = dt/2.0;
138                 for(int i=0;i<y.length-1;i++)
139                         y[i+1] = y[i] + dt*func.map(y[i]+dt2*func.map(y[i]));
140                 return y;
141         }
142         /**
143         * Uses the 2nd order Runge-Kutta method to solve an ODE.
144         * @param y an array to be filled with y values, set y[0] to initial condition.
145         * @param func dy/dt as a function of y and t.
146         * @param t0 initial time.
147         * @param dt step size.
148         * @return y.
149         */
150         public static double[] rungeKutta2(final double y[],final RealFunction2D func,final double t0, final double dt) {
151                 final double dt2 = dt/2.0;
152                 double t = t0;
153                 for(int i=0;i<y.length-1;i++) {
154                         y[i+1] = y[i] + dt*func.map(y[i]+dt2*func.map(y[i], t), t+dt2);
155                         t += dt;
156                 }
157                 return y;
158         }
159         /**
160         * Uses the 4th order Runge-Kutta method to solve an ODE.
161         * @param y an array to be filled with y values, set y[0] to initial condition.
162         * @param func dy/dt as a function of y.
163         * @param dt step size.
164         * @return y.
165         */
166         public static double[] rungeKutta4(final double y[],final Mapping func,final double dt) {
167                 double k1,k2,k3,k4;
168                 for(int i=0;i<y.length-1;i++) {
169                         k1 = dt*func.map(y[i]);
170                         k2 = dt*func.map(y[i]+k1/2.0);
171                         k3 = dt*func.map(y[i]+k2/2.0);
172                         k4 = dt*func.map(y[i]+k3);
173                         y[i+1] = y[i] + (k1+k4)/6.0 + (k2+k3)/3.0;
174                 }
175                 return y;
176         }
177         /**
178         * Uses the 4th order Runge-Kutta method to solve an ODE.
179         * @param y an array to be filled with y values, set y[0] to initial condition.
180         * @param func dy/dt as a function of y and t.
181         * @param dt step size.
182         * @return y.
183         */
184         public static double[] rungeKutta4(final double y[],final RealFunction2D func,final double t0, final double dt) {
185                 final double dt2 = dt/2.0;
186                 double t = t0;
187                 double k1,k2,k3,k4;
188                 for(int i=0;i<y.length-1;i++) {
189                         k1 = dt*func.map(y[i], t);
190                         k2 = dt*func.map(y[i]+k1/2.0, t+dt2);
191                         k3 = dt*func.map(y[i]+k2/2.0, t+dt2);
192                         k4 = dt*func.map(y[i]+k3, t+dt);
193                         y[i+1] = y[i] + (k1+k4)/6.0 + (k2+k3)/3.0;
194                         t += dt;
195                 }
196                 return y;
197         }
198         /**
199         * Numerical integration using the trapezium rule.
200         * @param N the number of strips to use.
201         * @param func a function.
202         * @param a the first ordinate.
203         * @param b the last ordinate.
204         */
205         public static double trapezium(final int N,final Mapping func,final double a,final double b) {
206                 double A=0.0,x=a,h=(b-a)/N;
207                 for(int i=0;i<N;i++) {
208                         A+=func.map(x)+func.map(x+h);
209                         x+=h;
210                 }
211                 return A*h/2.0;
212         }
213         /**
214         * Numerical integration using Simpson's rule.
215         * @param N the number of strip pairs to use.
216         * @param func a function.
217         * @param a the first ordinate.
218         * @param b the last ordinate.
219         */
220         public static double simpson(final int N,final Mapping func,final double a,final double b) {
221                 double Ao=0.0,Ae=0.0,x=a;
222                 final double h=(b-a)/(2*N);
223                 for(int i=0;i<N-1;i++) {
224                         Ao+=func.map(x+h);
225                         Ae+=func.map(x+2*h);
226                         x+=2.0*h;
227                 }
228                 Ao+=func.map(x+h);
229                 return h/3.0*(func.map(a)+4.0*Ao+2.0*Ae+func.map(b));
230         }
231         /**
232         * Numerical integration using the Richardson extrapolation.
233         * @param N the number of strip pairs to use (lower value).
234         * @param func a function.
235         * @param a the first ordinate.
236         * @param b the last ordinate.
237         */
238         public static double richardson(final int N,final Mapping func,final double a,final double b) {
239                 double Aa,Aao=0.0,Aae=0.0,Ab,Abo=0.0,Abe=0.0,x=a;
240                 final double ha=(b-a)/(2*N);
241                 final double hb=ha/2.0;
242                 for(int i=0;i<N-1;i++) {
243                         Aao+=func.map(x+ha);
244                         Aae+=func.map(x+2.0*ha);
245                         Abo+=func.map(x+hb);
246                         Abe+=func.map(x+2*hb);
247                         Abo+=func.map(x+3*hb);
248                         Abe+=func.map(x+4*hb);
249                         x+=2.0*ha;
250                 }
251                 Aao+=func.map(x+ha);
252                 Abo+=func.map(x+hb);
253                 Abe+=func.map(x+2.0*hb);
254                 Abo+=func.map(x+3.0*hb);
255                 Aa=ha/3.0*(func.map(a)+4.0*Aao+2.0*Aae+func.map(b));
256                 Ab=hb/3.0*(func.map(a)+4.0*Abo+2.0*Abe+func.map(b));
257                 return (16.0*Ab-Aa)/15.0;
258         }
259         /**
260         * Numerical integration using the Gaussian integration formula (4 points).
261         * @param N the number of strips to use.
262         * @param func a function.
263         * @param a the first ordinate.
264         * @param b the last ordinate.
265         */
266         public static double gaussian4(final int N,final Mapping func,double a,final double b) {
267                 int n,i;
268                 double A=0.0;
269                 final double h=(b-a)/N;
270                 final double h2=h/2.0;
271                 final double zeros[]=new double[4];
272                 final double coeffs[]=new double[4];
273                 zeros[2]=0.339981043584856264802665759103;
274                 zeros[3]=0.861136311594052575223946488893;
275                 zeros[0]=-zeros[3];
276                 zeros[1]=-zeros[2];
277                 coeffs[0]=coeffs[3]=0.347854845137453857373063949222;
278                 coeffs[1]=coeffs[2]=0.652145154862546142626936050778;
279                 for(n=0;n<N;n++) {
280                         for(i=0;i<zeros.length;i++)
281                                 A+=coeffs[i]*func.map(a+(zeros[i]+1)*h2);
282                         a+=h;
283                 }
284                 return A*h2;
285         }
286         /**
287         * Numerical integration using the Gaussian integration formula (8 points).
288         * @param N the number of strips to use.
289         * @param func a function.
290         * @param a the first ordinate.
291         * @param b the last ordinate.
292         */
293         public static double gaussian8(final int N,final Mapping func,double a,final double b) {
294                 int n,i;
295                 double A=0.0;
296                 final double h=(b-a)/N;
297                 final double h2=h/2.0;
298                 final double zeros[]=new double[8];
299                 final double coeffs[]=new double[8];
300                 zeros[4]=0.183434642495649804939476142360;
301                 zeros[5]=0.525532409916328985817739049189;
302                 zeros[6]=0.796666477413626739591553936476;
303                 zeros[7]=0.960289856497536231683560868569;
304                 zeros[0]=-zeros[7];
305                 zeros[1]=-zeros[6];
306                 zeros[2]=-zeros[5];
307                 zeros[3]=-zeros[4];
308                 coeffs[0]=coeffs[7]=0.101228536290376259152531354310;
309                 coeffs[1]=coeffs[6]=0.222381034453374470544355994426;
310                 coeffs[2]=coeffs[5]=0.313706645877887287337962201987;
311                 coeffs[3]=coeffs[4]=0.362683783378361982965150449277;
312                 for(n=0;n<N;n++) {
313                         for(i=0;i<zeros.length;i++)
314                                 A+=coeffs[i]*func.map(a+(zeros[i]+1)*h2);
315                         a+=h;
316                 }
317                 return A*h2;
318         }
319         /**
320         * Numerical differentiation.
321         * @param N the number of points to use.
322         * @param func a function.
323         * @param a the first ordinate.
324         * @param b the last ordinate.
325         */
326         public static double[] differentiate(final int N,final Mapping func,final double a,final double b) {
327                 final double diff[]=new double[N];
328                 double x=a;
329                 final double dx=(b-a)/N;
330                 final double dx2=dx/2.0;
331                 for(int i=0;i<N;i++) {
332                         diff[i]=(func.map(x+dx2)-func.map(x-dx2))/dx;
333                         x+=dx;
334                 }
335                 return diff;
336         }
337         /**
338         * Numerical differentiation in multiple dimensions.
339         * @param func a function.
340         * @param x coordinates at which to differentiate about.
341         * @param dx step size.
342         * @return an array M<sub>ij</sub>=df<sup>i</sup>/dx<sub>j</sub>.
343         */
344         public static double[][] differentiate(final MappingND func,final double x[],final double dx[]) {
345                 final double xplus[]=new double[x.length];
346                 final double xminus[]=new double[x.length];
347                 System.arraycopy(x,0,xplus,0,x.length);
348                 System.arraycopy(x,0,xminus,0,x.length);
349                 xplus[0]+=dx[0];
350                 xminus[0]-=dx[0];
351                 double funcdiff[]=ArrayMath.scalarMultiply(0.5/dx[0],ArrayMath.subtract(func.map(xplus),func.map(xminus)));
352                 final double diff[][]=new double[funcdiff.length][x.length];
353                 for(int i=0;i<funcdiff.length;i++)
354                         diff[i][0]=funcdiff[i];
355                 for(int i,j=1;j<x.length;j++) {
356                         System.arraycopy(x,0,xplus,0,x.length);
357                         System.arraycopy(x,0,xminus,0,x.length);
358                         xplus[j]+=dx[j];
359                         xminus[j]-=dx[j];
360                         funcdiff=ArrayMath.scalarMultiply(0.5/dx[j],ArrayMath.subtract(func.map(xplus),func.map(xminus)));
361                         for(i=0;i<funcdiff.length;i++)
362                                 diff[i][j]=funcdiff[i];
363                 }
364                 return diff;
365         }
366         /**
367         * The Metropolis algorithm.
368         * @param list an array to be filled with values distributed according to func, set list[0] to initial value.
369         * @param func distribution function.
370         * @param dx step size.
371         * @return list.
372         */
373         public static double[] metropolis(final double list[],final Mapping func,final double dx) {
374                 for(int i=0;i<list.length-1;i++) {
375                         list[i+1]=list[i]+dx*(2.0*Math.random()-1.0);
376                         if(func.map(list[i+1])/func.map(list[i]) < Math.random())
377                                 list[i+1]=list[i];
378                 }
379                 return list;
380         }
381 }
382 
383
A to Z: JavaDoc & Examples Daily Java News & Articles Open Source Projects Open Source Codes Free Computer Books Remove Frame
Popular Tags