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