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
|
/* Wronksian test for Bessel functions. */
long double jnl (), ynl (), floorl ();
#define PI 3.14159265358979323846L
long double y, Jn, Jnp1, Jmn, Jmnp1, Yn, Ynp1;
long double w1, w2, err1, max1, err2, max2;
void wronk ();
main ()
{
long double x, delta;
int n, i, j;
max1 = 0.0L;
max2 = 0.0L;
delta = 0.6 / PI;
for (n = -30; n <= 30; n++)
{
x = -30.0;
while (x < 30.0)
{
wronk (n, x);
x += delta;
}
delta += .00123456;
}
}
void
wronk (n, x)
int n;
long double x;
{
Jnp1 = jnl (n + 1, x);
Jmn = jnl (-n, x);
Jn = jnl (n, x);
Jmnp1 = jnl (-(n + 1), x);
/* This should be trivially zero. */
err1 = Jnp1 * Jmn + Jn * Jmnp1;
if (err1 < 0.0)
err1 = -err1;
if (err1 > max1)
{
max1 = err1;
printf ("1 %3d %.5Le %.3Le\n", n, x, max1);
}
if (x < 0.0)
{
x = -x;
Jn = jnl (n, x);
Jnp1 = jnl (n + 1, x);
}
Yn = ynl (n, x);
Ynp1 = ynl (n + 1, x);
/* The Wronksian. */
w1 = Jnp1 * Yn - Jn * Ynp1;
/* What the Wronksian should be. */
w2 = 2.0 / (PI * x);
err2 = w1 - w2;
if (err2 > max2)
{
max2 = err2;
printf ("2 %3d %.5Le %.3Le\n", n, x, max2);
}
}
|
