Added 3.4.2005: I've provided some hints about computing J1(x).
Tentatively Due before the beginning of class Thursday, March 9.
This assignment involves two Bessel functions. Bessel functions arise in studying physical systems best described in cylindrical coordinates, such as the flow of heat, current, or perhaps neutron flux, in a solid with a circular cross section. I first encountered them in my first work in computing, which was at a General Electric plant that produced plutonium in nuclear reactors in which the fuel was (in effect) a long rod of uranium.
We will not look at the differential equation from which Bessel functions arise, but only at their Taylor series representation. If you have not yet taken Math 203, then this will be your first introduction to Taylor series. I will give brief presentations of everything you need to know to do this assignment.
There are a great many Bessel functions. We will consider only two, J0 and J1, which are something like a damped cosine and a damped sine respectively. Here are the series, along with the basic transformations we make to speed computation. (If you can't read my handwriting, Google! Or see scan following.)
We will do this assignment in two stages:
To know whether we are getting the right answers, we need a table of values, and we have one. (This also gives code for doing the computation in Basic. The approach is the "obvious," "straightforward" way--which is unspeakably slow. We'll see why.)
I don't expect you to rediscover everything about this, in two weeks. Here's a starter kit of the computation and the testing for correctness. You will get more starter kit stuff later.
// DDM 2.19.2006 HW3: Bessel function experimentation
// Watch out! I used J0 as both the name of an array and the
// name of a method!!! Don't think I would ordinarily do this,
// but maybe once to make a point that the compiler has no problem.
package hw3;
public class Bessel {
// Array J0 of known correct values, for x = 0, 0.1, 0.2, ...
double J0[] = { 1, 0.997502, 0.990025, 0.977626, 0.960398, 0.93847,
0.912005, 0.881201, 0.846287, 0.807524, 0.765198, 0.719622,
0.671133, 0.620086 };
// Constructor
Bessel() {
for (int i = 0; i < J0.length; i++)
System.out.println(i / 10.0 + " " + J0[i] + " " + J0(i / 10.0));
}
// Method J0 for computing J0 from
// Taylor series, using Horner's method
// Starting value of n to get required accuracy
// depends on the size of the argument. This is part
// of the assignment.
double J0(double x) {
double sum = 1.0;
double y = x / 2.0;
double z = y * y;
for (int n = 6; n > 0; n--)
sum = 1 - (z / (n * n)) * sum;
return sum;
}
public static void main(String[] fred) {
new Bessel();
}
}
Added 3.4.2005: Here are some hints on writing the program for J1(x).
Let z = x/2, let ^ stand for exponentiation, and write:
J1(x) = z/(1) - z^3/(1*1*2) + z^5/(1*1*2*2*3) - z^7(1*1*2*2*3*3*4) + - ...
Stop with terms through z^7 and factor out z/(1):
J1(x) ~= (z/(1))*(1 - z^2/(1*2) + z^4/(1*2*2*3) - z^6/(1*2*2*3*3*4))
where ~= means "approximately equals."
Now look at the long expression in parentheses. Starting with the z^2/(1*2)
term, factor that out. Do the a similar kind of thing again, and you've
got the Horner's Rule version of the truncated series.
Now look at the program above and figure out how you need to initialize sum,
and what you need to do to the for-loop, to get J1.
JUnit/Eclipse is the Java name for unit testing, a methodology that says to write the test routine for a program first, and then write a program that does what the test routine expects. We don't do exactly that, but we learn the techniques.
A quick example will provide an overview of the subject. We have a class named Calculator with a method Sum that computes the sum of two doubles:
public class Calculator {
double Sum(double a, double b) {
return a + b;
}
}
We want to convince ourselves that the function works as intended. We write a class that extends a class named TestCase from a JUnit package:
import junit.framework.TestCase;
public class TestCalculator extends TestCase {
public void testAdd(){
Calculator calculator = new Calculator();
double expected = 12.0;
double computed = calculator.Sum(10.5, 1.5);
assertEquals(expected, computed);
}
}
This class tests a class named Calculator; by convention it is named TestCalculator. It contains a method that must be named testXXX, where XXX is anything we wish, to indicate what is being tested. We instantiate the Calculator class. Now we want to see whether a calculated sum is what we expect. We set a variable named expected to some value, then use the Sum method to compute a value that should be the same. (But watch out for comparing doubles for exact equality; we'll be back to this issue.)
Now comes the heart of the thing. assertEquals is a method from the class we are extending, TestCase. There are many overloaded versions of this; we pick one of the simplest, which compares two doubles for equality. This can be read as saying, "I assert that these two values are--or anyway should be--equal. Are they?"
Now we run . . . what? There is no main method anywhere here. New plan. When we click on Run As in Eclipse, we are offered the choice JUnit Test. We choose that, hold our breath, and see what happens. If things are working as they should (don't they always???) we see a beautiful green bar in a new tab at the left, called JUnit. If not, there is a red bar.
To get back to the Package Explorer view, just click on the small number (2 in this screenshot) to the right of the JUnit tab, and pull down to Package Explorer.
================
OK, I'm out of time for right now. Here's a class to compute the Bessel function J0, with two arguments: a value for x, and a value for the highest term number. This will let us experiment with how many terms we need to get desired accuracy.
public class BesselCalc {
double J0(double x, int highPower) {
double sum = 1.0;
double y = x / 2.0;
double z = y * y;
for (int n = highPower; n > 0; n--)
sum = 1 - (z / (n * n)) * sum;
return sum;
}
}
Now here is the JUnit class to do one test:
import junit.framework.TestCase;
public class TestBesselCalc extends TestCase {
public void testJ0Calc() {
BesselCalc besselCalc = new BesselCalc();
// Test J0(4.0)
double expected = -0.39715;
double computed = besselCalc.J0(4.0, 7);
assertEquals(expected, computed, 0.0001);
}
}
This isn't so different from the simple Sum test. The expected value comes from the table, and we call J0 for x = 4.0 and the highest value being 7, meaning terms through (x/2) raised to the 14th power.
The main difference is that we use a different overloaded assertEquals; this one takes two doubles to be compared, and a tolerance. The two values are considered equal if the absolute value of their difference is less than the tolerance.
Sketch of what's to come: You write a loop in TestBesselCalc that runs through values of x from 0.0 to 4.0, checking against table values.
When that all works correctly, you do the same thing for J1, using a Horner's rule form that you devise. Test it.
Now you can back off, and see how few terms you can get by with, when the argument is less than 2.0 and the tolerance is 0.0001, which is what our application requires.