KickJava   Java API By Example, From Geeks To Geeks.

Java > Open Source Codes > Integrate


1 import EDU.oswego.cs.dl.util.concurrent.*;
2
3
4 /**
5  * Sample program using Guassian Quadrature for numerical integration.
6  * Inspired by a
7  * <A HREF="http://www.cs.uga.edu/~dkl/filaments/dist.html"> Filaments</A>
8  * demo program.
9  *
10  */
11
12 public class Integrate {
13
14
15   public static void main(String JavaDoc[] args) {
16     int procs;
17     double start;
18     double end;
19     int exp = 5;
20
21     try {
22       procs = Integer.parseInt(args[0]);
23       start = new Double JavaDoc(args[1]).doubleValue();
24       end = new Double JavaDoc(args[2]).doubleValue();
25       if (args.length > 3) exp = Integer.parseInt(args[3]);
26     }
27     catch (Exception JavaDoc e) {
28       System.out.println("Usage: java Integrate <threads> <lower bound> <upper bound> <exponent>\n (for example 2 1 48 5).");
29       return;
30     }
31
32     System.out.println("Integrating from " + start + " to " + end + " exponent: " + exp);
33
34     Function f = new SampleFunction(exp);
35     FJTaskRunnerGroup group = new FJTaskRunnerGroup(procs);
36     Integrator integrator = new Integrator(f, 0.001, group);
37     double result = integrator.integral(start, end);
38     
39     System.out.println("Answer = " + result);
40     group.stats();
41
42   }
43
44   /*
45     This is all set up as if it were part of a more serious
46     framework, but is for now just a demo, with all
47     classes declared as static within Integrate
48   */
49
50   /** A function to be integrated **/
51   static interface Function {
52     double compute(double x);
53   }
54
55   /**
56    * Sample from filaments demo.
57    * Computes (2*n-1)*(x^(2*n-1)) for all odd values
58    **/
59   static class SampleFunction implements Function {
60     final int n;
61     SampleFunction(int n) { this.n = n; }
62
63     public double compute(double x) {
64       double power = x;
65       double xsq = x * x;
66       double val = power;
67       double di = 1.0;
68       for (int i = n - 1; i > 0; --i) {
69         di += 2.0;
70         power *= xsq;
71         val += di * power;
72       }
73       return val;
74     }
75   }
76
77
78   static class Integrator {
79     final Function f; // The function to integrate
80 final double errorTolerance;
81     final FJTaskRunnerGroup group;
82
83     Integrator(Function f, double errorTolerance, FJTaskRunnerGroup group) {
84       this.f = f;
85       this.errorTolerance = errorTolerance;
86       this.group = group;
87     }
88
89     double integral(double lowerBound, double upperBound) {
90       double f_lower = f.compute(lowerBound);
91       double f_upper = f.compute(upperBound);
92       double initialArea = 0.5 * (upperBound-lowerBound) * (f_upper + f_lower);
93       Quad q = new Quad(lowerBound, upperBound,
94                         f_lower, f_upper,
95                         initialArea);
96       try {
97         group.invoke(q);
98         return q.area;
99       }
100       catch(InterruptedException JavaDoc ex) {
101         Thread.currentThread().interrupt();
102         throw new Error JavaDoc("Interrupted during computation");
103       }
104     }
105
106     
107     /**
108      * FJTask to recursively perform the quadrature.
109      * Algorithm:
110      * Compute the area from lower bound to the center point of interval,
111      * and from the center point to the upper bound. If this
112      * differs from the value from lower to upper by more than
113      * the error tolerance, recurse on each half.
114      **/
115     class Quad extends FJTask {
116       final double left; // lower bound
117 final double right; // upper bound
118 final double f_left; // value of the function evaluated at left
119 final double f_right; // value of the function evaluated at right
120
121       // Area initialized with original estimate from left to right.
122 // It is replaced with refined value.
123 volatile double area;
124     
125       Quad(double left, double right,
126            double f_left, double f_right,
127            double area) {
128         this.left = left;
129         this.right = right;
130         this.f_left = f_left;
131         this.f_right = f_right;
132         this.area = area;
133       }
134       
135       public void run() {
136         double center = 0.5 * (left + right);
137         double f_center = f.compute(center);
138         
139         double leftArea = 0.5 * (center - left) * (f_left + f_center);
140         double rightArea = 0.5 * (right - center) * (f_center + f_right);
141         double sum = leftArea + rightArea;
142         
143         double diff = sum - area;
144         if (diff < 0) diff = -diff;
145         
146         if (diff >= errorTolerance) {
147           Quad q1 = new Quad(left, center, f_left, f_center, leftArea);
148           Quad q2 = new Quad(center, right, f_center, f_right, rightArea);
149           coInvoke(q1, q2);
150           sum = q1.area + q2.area;
151         }
152         
153         area = sum;
154       }
155     }
156   }
157
158 }
159
160   
161
Popular Tags