← Back to branch summary

~ubuntu-branches/debian/squeeze/maxima/squeeze

« back to all changes in this revision

Viewing changes to tests/rtest14.mac

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2006-10-18 14:52:42 UTC
  • mto: (1.1.5 upstream)
  • mto: This revision was merged to the branch mainline in revision 4.
  • Revision ID: james.westby@ubuntu.com-20061018145242-vzyrm5hmxr8kiosf
Import upstream version 5.10.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/rtest14.mac
1
1
jn(3,4);
2
 
0.1320341839226408 ;
 
2
jn(3,4); /*0.1320341839226408 ;*/
3
3
j0(1);
4
 
0.7651976865579665 ;
 
4
j0(1); /*0.7651976865579665 ;*/
5
5
bessel_j(0,1.0);
6
6
0.7651976865579665 ;
7
7
bessel(2,3);
 
8
bessel(2,3)$
 
9
bessel(2,3.0);
8
10
0.1289432494744021;
9
11
bessel_j(3,2.0);
10
12
0.1289432494744021;
91
93
true$
92
94
(assume(4*p+a>0),true);
93
95
true$
 
96
besselexpand:false;
 
97
false$
94
98
 
95
99
specint(t^(1/2)*%e^(-a*t/4)*%e^(-p*t),t);
96
100
sqrt(%pi)/(2*(p+a/4)^(3/2));
97
101
 
98
 
specint(t^(3/4)*%e^(-t^2/2/b)*%e^(-p*t),t);
99
 
3*gamma(3/4)*b^(7/8)*%e^(b*p^2/4)*(8*sqrt(%pi)*%m[-5/8,-1/4](b*p^2/2)/(3*2^ 
100
 
 (5/8)*gamma(3/8)*b^(1/4)*sqrt(p))-2^(3/8)*sqrt(%pi)*%m[-5/8,1/4](b*p^2/2)/ 
101
 
 (gamma(7/8)*b^(1/4)*sqrt(p)))/4 $
102
 
 
103
 
specint(t^(-1/2)*%e^(-2*a^(1/2)*t^(1/2))*%e^(-p*t),t);
104
 
sqrt(2)*%e^(a/(2*p))*(sqrt(%pi)*%e^(a/(2*p))/sqrt(2)-sqrt(2)*sqrt(%pi)*erf 
105
 
 (sqrt(a)/sqrt(p))*%e^(a/(2*p))/2)/sqrt(p) $
106
 
 
 
102
prefer_whittaker:true;
 
103
true$
 
104
 
 
105
/*
 
106
 * Reference:  Table of Integral Transforms.
 
107
 */
 
108
 
 
109
/*
 
110
 * f24p146:
 
111
 *
 
112
 * t^(v-1)*exp(-t^2/(8*a))
 
113
 *     -> gamma(v)*2^v*a^(v/2)*exp(a*p^2)*D[-v](2*p*sqrt(a))
 
114
 *
 
115
 * We have v = 7/4, a = b/4 so the result should be
 
116
 *
 
117
 *   gamma(7/4)*2^(7/4)*(b/4)^(7/8)*exp(b*p^2/4)*D[-7/4](p*sqrt(b))
 
118
 *     = 3/4*gamma(3/4)*b^(7/8)*exp(b*p^2/4)*D[-7/4](p*sqrt(b))
 
119
 *
 
120
 * But
 
121
 *
 
122
 *   D[v](z) = 2^(v/2+1/4)*z^(-1/2)*%w[v/2+1/4,1/4](z^2/2)
 
123
 *
 
124
 * and
 
125
 *
 
126
 *   %w[k,u](z) = gamma(-2*u)/gamma(1/2-u-k)*%m[k,u](z)
 
127
 *                 + gamma(2*u)/gamma(1/2+u-k)*%m[k,-u](z)
 
128
 *
 
129
 * Thus,
 
130
 *
 
131
 *   D[-7/4](p*sqrt(b)) = %w[-5/8,1/4](b*p^2/2)/(2^(5/8)*b^(1/4)*sqrt(p))
 
132
 *
 
133
 * and
 
134
 *
 
135
 *   %w[-5/8,1/4](b*p^2) 
 
136
 *     = 8/3/gamma(3/8)*%m[-5/8,-1/4](b*p^2/2)
 
137
 *         - 2*sqrt(%pi)/gamma(7/8)*%m[-5/8,1/4](b*p^2/2)
 
138
 *
 
139
 * And finally the transform is
 
140
 *
 
141
 * 3*gamma(3/4)*b^(5/8)*exp(b*p^2/4)/(4*2^(5/8)*sqrt(p))
 
142
 *   *(8*sqrt(%pi)/3/gamma(3/8)*%m[-5/8,-1/4](b*p^2/2)
 
143
 *      - 2*sqrt(%pi)/gamma(7/8)*%m[-5/8,1/4](b*p^2/2))
 
144
 */
 
145
(assume(b>0),true);
 
146
true$
 
147
 
 
148
ratsimp(specint(t^(3/4)*%e^(-t^2/2/b)*%e^(-p*t),t));
 
149
-sqrt(%pi)*b^(5/8)
 
150
          *(3*gamma(3/8)*gamma(3/4)*%e^(b*p^2/4)*%m[-5/8,1/4](b*p^2/2)
 
151
           -4*gamma(3/4)*gamma(7/8)*%e^(b*p^2/4)
 
152
             *%m[-5/8,-1/4](b*p^2/2))
 
153
 /(2*2^(5/8)*gamma(3/8)*gamma(7/8)*sqrt(p))$
 
154
 
 
155
/*
 
156
 * Sec. 4.5, formula (33):
 
157
 *
 
158
 * t^(-1/2)*exp(-2*sqrt(a)*sqrt(t)) ->
 
159
 *    sqrt(%pi)/sqrt(p)*exp(a/p)*erfc(sqrt(a)/sqrt(p))
 
160
 */
 
161
ratsimp(specint(t^(-1/2)*%e^(-2*a^(1/2)*t^(1/2))*%e^(-p*t),t));
 
162
-sqrt(%pi)*(erf(sqrt(a)/sqrt(p))-1)*%e^(a/p)/sqrt(p)$
 
163
 
 
164
/*
 
165
 * The Laplace transform of sin(a*t)*cosh(b*t^2) can be derived by
 
166
 * expressing sin and cosh in terms of exponential functions.  We end
 
167
 * up with terms of the form:
 
168
 *
 
169
 *   exp(+/-%i*a*t)*exp(+/-b*t^2)
 
170
 *
 
171
 * All of these can be computed using formula 24, p. 146 of Tables of
 
172
 * Integral Transforms, which handle functions of the form
 
173
 * t^(v-1)*exp(-t^2/8/a).
 
174
 *
 
175
 * But, we have terms of the form exp(b*t^2-p*t+%i*a*t).  I don't
 
176
 * think this converges, so the Laplace transform doesn't exist if b >
 
177
 * 0.
 
178
 * 
 
179
 */
 
180
/*
107
181
radcan(specint(sin(a*t)*cosh(b*t^2)*%e^(-p*t),t));
108
182
-%e^-((p^2+2*%i*a*p+a^2)/(4*b))*(sqrt(%pi)*%e^((2*%i*a*p+a^2)/(2*b))*erf((%i 
109
183
 *p+a)/(2*sqrt(b)))-sqrt(%pi)*%e^(a^2/(2*b))*erf((%i*p-a)/(2*sqrt(b)))+sqrt 
111
185
 ^(p^2/(2*b))*erf((p-%i*a)/(2*sqrt(b)))+(sqrt(%pi)*%i-sqrt(%pi)*%i*%e^(%i*a 
112
186
 *p/b))*%e^(p^2/(2*b))-sqrt(%pi)*%e^((2*%i*a*p+a^2)/(2*b))+sqrt(%pi)*%e^(a^2/ 
113
187
 (2*b)))/(8*sqrt(b)) $
114
 
 
115
 
 
116
 
specint(t^(1/2)*%j[1](2*a^(1/2)*t^(1/2))*%e^(-p*t),t);
 
188
*/
 
189
 
 
190
/*
 
191
 * Sec 4.14, formula (27):
 
192
 *
 
193
 * t^(1/2)*bessel_j(1,2*a^(1/2)*t^(1/2)) ->
 
194
 *    sqrt(a)/p^2*exp(-a/p)
 
195
 */
 
196
 
 
197
specint(t^(1/2)*bessel_j(1,2*a^(1/2)*t^(1/2))*%e^(-p*t),t);
117
198
sqrt(a)*%e^-(a/p)/p^2$
118
199
 
119
 
specint(t^2*%j[1](a*t)*%e^(-p*t),t);
120
 
3*a/((a^2/p^2+1)^(5/2)*p^4) ;
121
 
specint(t*hstruve[1](t)*%e^(-p*t),t);
122
 
16*(2*p*sqrt(p^2+1)+2*p^2)/(9*%pi^(3/2)*p^4*(2*p*sqrt(p^2+1)+2*p^2+1)) $
123
 
radcan(specint(t^(3/2)*hstruve[1](t^(1/2))*%e^(-p*t),t));
124
 
%e^-(1/(4*p))*((140*p^(3/2)+10*sqrt(p))*%e^(1/(4*p))+140*sqrt(%pi)*%i*erf(%i/ 
125
 
 (2*sqrt(p)))*p^2+20*sqrt(%pi)*%i*erf(%i/(2*sqrt(p)))*p+5*sqrt(%pi)*%i*erf 
126
 
 (%i/(2*sqrt(p))))/(8*%pi*p^4) $
127
 
 
 
200
/*
 
201
 * Sec 4.14, formula (3):
 
202
 *
 
203
 * t^2*bessel_j(v,a*t) ->
 
204
 *    ((v^2-1)/r^3 + 3*p*(p+v*r)/r^5)*(a/R)^v
 
205
 *
 
206
 * where r = sqrt(p^2+a^2), R = p + r;
 
207
 *
 
208
 * (Maxima can't currently compute this transform for general v due to a bug
 
209
 * in hyp.lisp.)
 
210
 */
 
211
factor(ratsimp(specint(t^2*bessel_j(1,a*t)*%e^(-p*t),t)));
 
212
3*a*p/(p^2+a^2)^(5/2) $
 
213
 
 
214
(/* This is the Laplace transform of the Struve H_1 function, see
 
215
  http://dlmf.nist.gov/Draft/ST/about_ST.8.13.html */    
 
216
 2/(%pi*p)-2*p*log(p/(sqrt(p^2+1)-1))/(%pi*sqrt(p^2+1)),
 
217
 /* And this should be the same as the specint of the next test below */
 
218
 -diff(%%,p),
 
219
 ev(fullratsimp(%%),logexpand:all));
 
220
-(sqrt(p^2+1)*(2*p^2*log(sqrt(p^2+1)-1)-2*p^2*log(p))-2*p^2-2)
 
221
 /(%pi*p^6+2*%pi*p^4+%pi*p^2)$
 
222
 
 
223
(ev(fullratsimp(specint(t*hstruve[1](t)*%e^(-p*t),t)),logexpand:all),
 
224
 ratsimp(%%/%));
 
225
1$
 
226
 
 
227
/*
 
228
 *
 
229
 * From the comments for hstf in hypgeo.lisp:
 
230
 *
 
231
 * hstruve[1](t) = 2/sqrt(%pi)*(t/2)^2/gamma(1+3/2)*%f[1,2]([1],[3/2,5/2],-t^2/4)
 
232
 *
 
233
 * So
 
234
 *
 
235
 * hstruve[1](sqrt(t)) = 2/(3*%pi)*t*%f[1,2]([1],[3/2,5/2],-t/4)
 
236
 *
 
237
 * and the integrand is
 
238
 *
 
239
 * 2/(3*%pi)*t^(5/2)*exp(-p*t)*%f[1,2]([1],[3/2,5/2],-t/4).
 
240
 *
 
241
 * From the f19p220, the Laplace transform of this, with s = 7/2,
 
242
 * c=-1/4, k = 1, is
 
243
 *
 
244
 * 2/(3*%pi)*gamma(7/2)/p^(7/2)*%f[2,2]([1,7/2],[3/2,5/2],-1/4/p)
 
245
 *
 
246
 * From the derivation of SPLITPFQ, we can simplify this
 
247
 * hypergeometric function.
 
248
 *
 
249
 * %f[2,2]([1,7/2],[3/2,5/2],z) =
 
250
 *
 
251
 *      1
 
252
 *     sum z^k/poch(5/2,k)*binomial(1,k) *diff(%f[2,2]([1,5/2],[3/2,5/2],z,k)
 
253
 *     k=0
 
254
 * 
 
255
 * But %f[2,2]([1,5/2],[3/2,5/2],z) = %f[1,1]([1],[3/2],z) 
 
256
 * and Maxima knows how to compute this.
 
257
 */
 
258
ratsimp(specint(t^(3/2)*hstruve[1](t^(1/2))*%e^(-p*t),t));
 
259
-%e^-(1/(4*p))*(sqrt(%pi)*sqrt(p)
 
260
                                *(8*%i*erf(%i/(2*sqrt(p)))*p
 
261
                                 -%i*erf(%i/(2*sqrt(p))))
 
262
                      -2*p*%e^(1/(4*p)))
 
263
        /(8*sqrt(%pi)*p^(9/2)) $
 
264
 
 
265
/* Trivial result because %ibes is not bessel_i
 
266
 */ 
128
267
specint(t*%ibes[0](a*t/2)*%ibes[1](a*t/2)*%e^(-p*t),t);
129
268
%ibes[0](a*t/2)*%ibes[1](a*t/2)/p^2 $
130
269
 
131
 
specint(t^(3/4)*%j[1/2](t)*%j[1/4](t)*%e^(-p*t),t);
 
270
/*
 
271
 * t^(3/4)*bessel_j(1/2,t)*bessel_j(1/4,t)
 
272
 *
 
273
 * Luke gives
 
274
 *
 
275
 * bessel_j(u,t)*bessel_j(v,t)
 
276
 *    = (z/2)^(u+v)/gamma(u+1)/gamma(v+1)
 
277
 *        * %f[2,3]([(u+v+1)/2,(u+v+2)/2],[u+1,v+1,u+v+1],-z^2)
 
278
 *
 
279
 * So the integrand is
 
280
 *
 
281
 * 8/2^(3/4)/sqrt(%pi)/gamma(1/4)*t^(3/2)*%f[2,3]([7/8,11/8],[3/2,5/4,7/4],-t^2)
 
282
 *
 
283
 * f19p220 gives
 
284
 *
 
285
 * t^(3/2)*%f[2,3]([7/8,11/8],[3/2,5/4,7/4],-t^2)
 
286
 *
 
287
 *    -> gamma(5/2)*p^(-5/2)
 
288
 *         *%f[4,3]([7/8,11/8,5/2/2,(5/2+1)/2],[3/2,5/4,7/4],-4/p^2)
 
289
 *    =  gamma(5/2)*p^(-5/2)
 
290
 *         $%f[2,1]([7/8,11/8],[3/2],-4/p^2)
 
291
 *
 
292
 * And we know %f[2,1]([7/8,11/8],[3/2],z) is
 
293
 *
 
294
 * -2*(1/(sqrt(z)+1)^(3/4)-1/(1-sqrt(z))^(3/4))/(3*sqrt(z))
 
295
 *
 
296
 * Applying all of this gives the expected answer below.
 
297
 */
 
298
 
 
299
specint(t^(3/4)*bessel_j(1/2,t)*bessel_j(1/4,t)*%e^(-p*t),t);
132
300
2*%i*(1/(2*%i/p+1)^(3/4)-1/(1-2*%i/p)^(3/4))/(2^(3/4)*gamma(1/4)*p^(3/2)) $
133
301
 
134
 
specint(t^(5/2)*%y[1/2](t^(1/2))^2*%e^(-p*t),t);
 
302
/*
 
303
 * Not sure this is right.  We can convert bessel_y to bessel_j,
 
304
 * multiply them together and use the results for products of bessel_j
 
305
 * functions.
 
306
 *
 
307
 * bessel_y(1/2,sqrt(t)) = -bessel_j(-1/2,sqrt(t))
 
308
 *
 
309
 * And maxima should be able to compute the transform of 
 
310
 *
 
311
 * t^(5/2)*bessel_j(-1/2,sqrt(t))^2
 
312
 *
 
313
 * Or note that bessel_y(1/2,sqrt(t)) =
 
314
 * -sqrt(2/%pi)*cos(sqrt(t))/t^(1/4).  Then the integrand becomes
 
315
 *
 
316
 * 2/%pi*t^2*cos(sqrt(t))^2
 
317
 *
 
318
 * And maxima should know how to compute the transform of this.
 
319
 *
 
320
 * Unfortunately, the transforms of these two approaches don't agree.
 
321
 * Yuck!
 
322
 */
 
323
specint(t^(5/2)*bessel_y(1/2,t^(1/2))^2*%e^(-p*t),t);
135
324
-12*((-3*sqrt(%pi)*%i*erf(%i/sqrt(p))*p^(5/2)*%e^-(1/p)/8-sqrt(%pi)*%i*erf 
136
325
 (%i/sqrt(p))*p^(3/2)*%e^-(1/p)/2-sqrt(%pi)*%i*erf(%i/sqrt(p))*sqrt(p)*%e^- 
137
326
 (1/p)/2-3*p^2/4-p/2)/p^2+2*(sqrt(%pi)*%i*erf(%i/sqrt(p))*p^(3/2)*%e^-(1/p)/4 
138
327
 +sqrt(%pi)*%i*erf(%i/sqrt(p))*sqrt(p)*%e^-(1/p)/2+p/2)/p-sqrt(%pi)*%i*erf 
139
328
 (%i/sqrt(p))*sqrt(p)*%e^-(1/p)/2)/(%pi*p^4) $
140
329
 
141
 
specint(t^(1/2)*%j[1/2](t^(1/2))^2*%e^(-p*t),t);
 
330
/*
 
331
 * See formula (42), p. 187:
 
332
 *
 
333
 * t^(lam-1)*bessel_j(2*u,2*sqrt(a)*sqrt(t))*bessel_j(2*v,2*sqrt(a)*sqrt(t))
 
334
 *
 
335
 *    -> 2*gamma(lam+u+v)*a^(u+v)/gamma(2*u+1)/gamma(2*v+1)/p^(lam+u+v)
 
336
 *        *%f[3,3]([u+v+1/2,u+v+1,lam+u+v],[2*u+1,2*v+1,2*u+2*v+1],-4*a/p)
 
337
 *
 
338
 * with Re(lam + u + v) > 0.
 
339
 *
 
340
 * So, we have lam = 3/2, u=v=1/4, a = 1/4, we get
 
341
 *
 
342
 * 4/%pi/p^2*%f[3,3]([1,3/2,2],[3/2,3/2,2],-1/p)
 
343
 *  = 4/%pi/p^2*%f[1,1]([1],[3/2],-1/p)
 
344
 *
 
345
 * And %f[1,1]([1],[3/2],-1/p) is 
 
346
 *
 
347
 *     -sqrt(%pi)*%i*erf(%i*sqrt(1/p))*%e^-(1/p)/(2*sqrt(1/p))
 
348
 *
 
349
 * So, the final result is:
 
350
 *
 
351
 * -2*%i*erf(%i/sqrt(p))*%e^-(1/p)/(sqrt(%pi)*p^(3/2))
 
352
 *
 
353
 * But we also have 
 
354
 *
 
355
 * bessel_j(u,t)*bessel_j(v,t)
 
356
 *    = (z/2)^(u+v)/gamma(u+1)/gamma(v+1)
 
357
 *        * %f[2,3]([(u+v+1)/2,(u+v+2)/2],[u+1,v+1,u+v+1],-z^2)
 
358
 *
 
359
 * So bessel_j(1/2,sqrt(t))^2 is
 
360
 *
 
361
 *    2/%pi*%f[2,3]([1,3/2],[3/2,3/2,2],-t)*sqrt(t)
 
362
 * 
 
363
 * So the integrand is
 
364
 *
 
365
 *    2/%pi*t*%f[2,3]([1,3/2],[3/2,3/2,2],-t)
 
366
 *     = 2/%pi*t*%f[1,2]([1],[3/2,2],-t)
 
367
 *
 
368
 * f19p220 then gives us the desired transform:
 
369
 *
 
370
 *    t*%f[1,2]([1],[3/2,2],-t)
 
371
 *      -> gamma(2)*p^(-2)*%f[2,2]([1,2],[3/2,2],-1/p)
 
372
 *
 
373
 *      = p^(-2)*%f[1,1]([1],[3/2],-1/p)
 
374
 *
 
375
 * So the final answer is
 
376
 *
 
377
 *    -%i*erf(%i/sqrt(p))*%e^-(1/p)/(sqrt(%pi)*p^(3/2))
 
378
 *
 
379
 * Hmm.  This differs from formula 42 above.  I think there's a bug in
 
380
 * formula 42, and it should be divided by 2.
 
381
 *
 
382
 * If we use the expression for the product of Bessel functions and
 
383
 * f19p220, we can easily derive the result of formula 42, except,
 
384
 * we're missing the factor of 2.  So, I think formula 42 is wrong.
 
385
 */
 
386
 
 
387
specint(t^(1/2)*bessel_j(1/2,t^(1/2))^2*%e^(-p*t),t);
142
388
-%i*erf(%i/sqrt(p))*%e^-(1/p)/(sqrt(%pi)*p^(3/2)) $
143
 
specint(t^(1/2)*%ibes[1](t)*%e^(-p*t),t);
144
 
sqrt(%pi)*%ibes[1](t)/(2*p^(3/2)) $
145
 
 
 
389
 
 
390
/*
 
391
 * See formula (8), section 4.16 of Table of Integral Transforms:
 
392
 *
 
393
 * t^u*bessel_i(v,a*t) -> gamma(u+v+1)*s^(-u-1)*assoc_legendre_p(u,-v,p/s)
 
394
 *
 
395
 * where s = sqrt(p^2-a^2);
 
396
 */
 
397
factor(ratsimp(specint(t^(1/2)*bessel_i(1,t)*%e^(-p*t),t)));
 
398
3*sqrt(%pi)*assoc_legendre_p(1/2,-1,p/sqrt(p^2-1))/(4*(p^2-1)^(3/4))$
 
399
 
 
400
/*
 
401
 * %h[2/3,1](sqrt(t)) = bessel_j(2/3,sqrt(t))+%i*bessel_y(2/3,sqrt(t))
 
402
 *
 
403
 * Formula (34) below gives:
 
404
 *
 
405
 * t^(u-1/2)*bessel_j(2*v,2*sqrt(a)*sqrt(t)) ->
 
406
 *   gamma(u+v+1/2)/sqrt(a)/gamma(2*v+1)*p^(-u)*exp(-a/p/2)*%m[u,v](a/p)
 
407
 *
 
408
 * Formula (50) gives
 
409
 *
 
410
 * t^(u-1/2)*bessel_y(2*v,2*sqrt(a)*sqrt(t)) ->
 
411
 *   1/sqrt(a)*p^(-u)*exp(-a/p/2)*
 
412
 *     (tan((u-v)*%pi)*gamma(u+v-1/2)/gamma(2*v+1) * %m[u,v](a/p)
 
413
 *       - sec((u-v)*%pi)*%w[u,v](a/p))
 
414
 *
 
415
 * But A&S 13.1.34 says
 
416
 *
 
417
 * %w[k,u](z) = gamma(-2*u)/gamma(1/2-u-k)*%m[k,u](z) 
 
418
 *               + gamma(2*u)/gamma(1/2+u-k)*%m[k,-u](z)
 
419
 *
 
420
 */
146
421
specint( t*%h[2/3,1](t^(1/2))*%e^(-p*t),t);
147
422
-4*%i*gamma(1/3)*%m[-3/2,1/3](-1/(4*p))*%e^-(1/(8*p))/(3*(-1)^(5/6)*sqrt(3) 
148
423
 *gamma(2/3)*p^(3/2))+4*gamma(1/3)*%m[-3/2,1/3](-1/(4*p))*%e^-(1/(8*p))/(3*( 
149
424
 -1)^(5/6)*gamma(2/3)*p^(3/2))-8*%i*gamma(2/3)*%m[-3/2,-1/3](-1/(4*p))*%e^- 
150
425
 (1/(8*p))/(3*(-1)^(1/6)*sqrt(3)*gamma(1/3)*p^(3/2)) $
151
426
 
 
427
/*
 
428
 * H[3/4,2](t) = bessel_j(3/4,t)-%i*bessel_y(3/4,t)
 
429
 *
 
430
 * Sec 4.14, formula (9):
 
431
 *
 
432
 * t^u*bessel_j(v,a*t) -> gamma(u+v+1)*r^(-u-1)*assoc_legendre_p(u,-v,p/r)
 
433
 *
 
434
 * where r = sqrt(p^2+a^2)
 
435
 *
 
436
 * Sec 4.14, formula (48)
 
437
 *
 
438
 * t^u*bessel_y(v,a*t) 
 
439
 *    -> r^(-u-1)*(gamma(u+v+1)*cot(v*%pi)*assoc_legendre_p(u,-v,p/r)
 
440
 *                  -gamma(u-v+1)*csc(v*%pi)*assoc_legendre_p(u,v,p/r))
 
441
 *
 
442
 * So, 
 
443
 *
 
444
 * t^(1/2)*bessej_j(3/4,t) 
 
445
 *    -> gamma(9/4)*r^(-3/2)*assoc_legendre_p(1/2,-3/4,p/r)
 
446
 *    =  5*gamma(1/4)/16*r^(-3/2)*assoc_legendre_p(1/2,-3/4,p/r)
 
447
 *
 
448
 * t^(1/2)*bessel_y(3/4,t) 
 
449
 *    -> r^(-3/2)*(gamma(9/4)*cot(3/4*%pi)*assoc_legendre_p(1/2,-3/4,p/r)
 
450
 *                  -gamma(3/4)*csc(3/4*%pi)*assoc_legendre_p(1/2,3/4,p/r))
 
451
 *    =  r^(-3/2)*(-gamma(9/4)*assoc_legendre_p(1/2,-3/4,p/r)
 
452
 *                  -gamma(3/4)*sqrt(2)*assoc_legendre_p(1/2,3/4,p/r))
 
453
 */
152
454
specint( t^(1/2)*%h[3/4,2](t)*%e^(-p*t),t);
153
 
5*%i*gamma(1/4)*%p[-3/2,-3/4](1/sqrt(1/p^2+1))*(1/p^4-1)^(3/8)*p^(9/4)/(18*( 
154
 
 -1)^(1/4)*sqrt(2)*gamma(3/4)^2)+5*gamma(1/4)*%p[-3/2,-3/4](1/sqrt(1/p^2+1))* 
155
 
 (1/p^4-1)^(3/8)*p^(9/4)/(18*(-1)^(1/4)*sqrt(2)*gamma(3/4)^2)+4*%i*gamma(3/4) 
156
 
 *%p[-3/2,3/4](1/sqrt(1/p^2+1))*p^(3/4)/((-1)^(3/4)*gamma(1/4)^2*(1/p^4-1)^ 
157
 
 (3/8)) $
 
455
5*gamma(1/4)/16/(p^2+1)^(3/4)*assoc_legendre_p(1/2,-3/4,p/sqrt(p^2+1))
 
456
+sqrt(2)*%i*assoc_legendre_p(1/2,3/4,p/sqrt(p^2+1))*gamma(3/4)
 
457
       /(p^2+1)^(3/4)
 
458
       +5*%i*gamma(1/4)*assoc_legendre_p(1/2,-3/4,p/sqrt(p^2+1))
 
459
        /(16*(p^2+1)^(3/4))$
158
460
 
 
461
/*
 
462
 * %h[1/2,1](t) = bessel_j(1/2,t)+%i*bessel_y(1/2,t)
 
463
 *
 
464
 * So,
 
465
 *
 
466
 * t^(3/2)*bessel_j(1/2,t)
 
467
 *     -> gamma(3/2+1/2+1)*r^(-5/2)*assoc_legendre_p(3/2,-1/2,p/r)
 
468
 *     =  2*r^(-5/2)*assoc_legendre_p(3/2,-1/2,p/r)
 
469
 * t^(3/2)*bessel_y(1/2,t) 
 
470
 *     -> r^(-5/2)*(gamma(3/2+1/2+1)*cot(%pi/2)*assoc_legendre_p(3/2,-1/2,p/r)
 
471
 *                  -gamma(3/2-1/2+1)*csc(%pi/2)*assoc_legendre_p(3/2,1/2,p/r))
 
472
 *     =  -r^(-5/2)*assoc_legendre_p(3/2,1/2,p/r))
 
473
 *
 
474
 * assoc_legendre_p(3/2,+/-1/2,z) can be expressed in terms of 
 
475
 * hypergeometric functions (A&S 8.1.2).
 
476
 *
 
477
 * assoc_legendre_p(3/2,-1/2,z) 
 
478
 *   = 1/gamma(3/2)*((z+1)/(z-1))^(-1/4)*F(-3/2,5/2;3/2;(1-z)/2)
 
479
 *   = sqrt(2)*(z-1)^(1/4)*z*(z+1)^(1/4)/sqrt(%pi)
 
480
 *
 
481
 * assoc_legendre_p(3/2,1/2,z) 
 
482
 *   = 1/gamma(-1/2)*((z+1)/(z-1))^(1/4)*F(-3/2,5/2;-1/2;(1-z)/2)
 
483
 *   = sqrt(2)*z*(2*z^2-3)/(sqrt(%pi)*(z-1)^(1/4)*(z+1)^(5/4))
 
484
 *
 
485
 *
 
486
 * So the result should be
 
487
 * 
 
488
 * t^(3/2)*bessel_j(1/2,t)
 
489
 *    -> 4*p/(sqrt(2)*sqrt(%pi)*(p^2+1)^2)
 
490
 *
 
491
 * t^(3/2)*bessel_y(1/2,t)
 
492
 *    -> -sqrt(2)*(p-1)*(p+1)/(sqrt(%pi)*(p^2+1)^2)
 
493
 */
159
494
specint( t^(3/2)*%h[1/2,1](t)*%e^(-p*t),t);
160
 
4/(sqrt(2)*sqrt(%pi)*(1/p^2+1)^2*p^3)-sqrt(2)*%i*(1/(1/p^2+1)-1/((1/p^2+1)^2 
161
 
 *p^2))/(sqrt(%pi)*p^2) $
162
 
 
163
 
specint( t^(3/2)*%y[1](a*t)*%e^(-t),t);
 
495
4/(sqrt(2)*sqrt(%pi)*(1/p^2+1)^2*p^3)
 
496
         -sqrt(2)*%i*(1/(1/p^2+1)-2/((1/p^2+1)^2*p^2))/(sqrt(%pi)*p^2) $
 
497
 
 
498
 
 
499
/*
 
500
 * Formula 2, p 105:
 
501
 *
 
502
 * t^(u-1)*bessel_y(v,a*t)
 
503
 *    -> -2/%pi*gamma(u+v)*(a^2+p^2)^(-u/2)
 
504
 *         *assoc_legendre_q(u-1,-v,p/sqrt(a^2+p^2))
 
505
 *
 
506
 * for a > 0, Re u > |Re v|
 
507
 *
 
508
 * We have u = 5/2, v = 1, so the result is
 
509
 *
 
510
 *    -4/%pi/(p^2+a^2)*assoc_legendre_q(3/2,-1,p/sqrt(p^2+a^2))
 
511
 *
 
512
 *
 
513
 */
 
514
specint( t^(3/2)*bessel_y(1,a*t)*%e^(-t),t);
164
515
15*%i*(1/(a^2+1)-1)^(3/4)*(1/(sqrt(a^2+1)+1)^(3/2)+1/(1-sqrt(a^2+1))^(3/2))/ 
165
516
 (16*sqrt(a^2+1)) $
166
517
 
167
 
specint( t^2*%j[1](a*t)*%e^(-p*t),t);
168
 
3*a/((a^2/p^2+1)^(5/2)*p^4) $
169
 
 
170
 
specint(t^(1/2)*%j[1](2*a^(1/2)*t^(1/2))*%e^(-p*t),t);
171
 
sqrt(a)*%e^-(a/p)/p^2 $
172
 
 
173
 
specint( t^(3/2)*%m[1/2,1](t)*%e^(-p*t),t);
174
 
6*(1/(1-1/(p+1/2))+1/((p+1/2)*(1-1/(p+1/2))^2))/(p+1/2)^4 ;
 
518
 
 
519
/*
 
520
 * A&S 13.1.32:
 
521
 *
 
522
 *   %m[k,u](t) = exp(-t/2)*t^(u+1/2)*M(1/2+u-k,1+2*u,t)
 
523
 *
 
524
 * So 
 
525
 *
 
526
 *   %m[1/2,1](t) = exp(-t/2)*t^(3/2)*M(1,3,t)
 
527
 *
 
528
 * and the integrand is:
 
529
 *
 
530
 *   t^(3/2)*%m[1/2,1](t)*exp(-p*t)
 
531
 *      = t^3*M(1,3,t)*exp(-(p+1/2)*t)
 
532
 *      = t^3*M(1,3,t)*exp(-p'*t)
 
533
 *
 
534
 * f19p220 will give us the Laplace transform of t^3*M(1,3,t):
 
535
 *
 
536
 *   gamma(4)/p'^4*F(1,4;3;1/p')
 
537
 * 
 
538
 * But p' = p+1/2, so the final result is
 
539
 *
 
540
 *   32*(6*p-1)/(2*p-1)^2/(2*p+1)^3
 
541
 * 
 
542
 */
 
543
ratsimp(specint( t^(3/2)*%m[1/2,1](t)*%e^(-p*t),t));
 
544
32*(6*p-1)/(2*p-1)^2/(2*p+1)^3;
 
545
 
175
546
(assume(p>a),true);
176
547
true;
177
 
 
 
548
/*
 
549
 * exp(a*t)*t^2*erf(sqrt(t))*exp(-p*t)
 
550
 * 
 
551
 * A&S 7.1.21 gives erf(z) = 2/sqrt(%pi)*z*M(1/2,3/2,-z^2) so 
 
552
 * erf(sqrt(t)) = 2/sqrt(%pi)*sqrt(t)*M(1/2,3/2,-t)
 
553
 *
 
554
 * Therefore, the integrand, with p' = p-a, is
 
555
 *
 
556
 *   2/sqrt(%pi)*t^(5/2)*M(1/2,3/2,-t)*exp(-p'*t)
 
557
 *
 
558
 * Applying f19p220, the Laplace transform is
 
559
 *
 
560
 *   2/sqrt(%pi)*gamma(7/2)/p'^(7/2)*F(1/2,7/2;3/2;-1/p')
 
561
 *
 
562
 * Maxima can compute F(1/2,7/2;3/2;-1/p') and substituting p'=p-a
 
563
 * gives us the desired answer.
 
564
 *
 
565
 * 15*(1/sqrt(1/(p-a)+1)-2/(3*(p-a)*(1/(p-a)+1)^(3/2))
 
566
 *                      +1/(5*(p-a)^2*(1/(p-a)+1)^(5/2)))
 
567
 *  /(4*(p-a)^(7/2))
 
568
 */
178
569
specint(%e^(a*t)*t^2*erf(t^(1/2))*%e^(-p*t),t);
179
 
15*(1/sqrt(1/(p-a)+1)-1/((p-a)*(1/(p-a)+1)^(3/2))+3/(4*(p-a)^2*(1/(p-a)+1)^ 
180
 
 (5/2)))/(4*(p-a)^(7/2)) ;
181
 
 
182
 
 
183
 
 
184
 
 
 
570
15*(1/sqrt(1/(p-a)+1)-2/(3*(p-a)*(1/(p-a)+1)^(3/2))
 
571
                     +1/(5*(p-a)^2*(1/(p-a)+1)^(5/2)))
 
572
 /(4*(p-a)^(7/2)) $
 
573
 
 
574
/*
 
575
 * Laplace transforms from Tables of Integral Transforms
 
576
 */
 
577
 
 
578
/*
 
579
 * p 182, (1)
 
580
 *
 
581
 * bessel_j(v,a*t) -> r^(-1)*(a/R)^v
 
582
 *
 
583
 * where r = sqrt(p^2+a^2) and R = p + r
 
584
 */
 
585
(assume(v>0),true);
 
586
true$
 
587
 
 
588
radcan(specint(bessel_j(v,a*t)*exp(-p*t),t));
 
589
a^v/(sqrt(p^2+a^2)*(sqrt(p^2+a^2)+p)^v)$
 
590
 
 
591
/*
 
592
 * (5)
 
593
 * bessel_j(v,a*t)/t -> v^(-1)*(a/R)^v
 
594
 *
 
595
 * (Maxima doesn't recognize that gamma(v)/gamma(v+1) is 1/v.)
 
596
 */
 
597
radcan(specint(bessel_j(v,a*t)/t*exp(-p*t),t));
 
598
a^v*gamma(v)/((sqrt(p^2+a^2)+p)^v*gamma(v+1))$
 
599
 
 
600
/*
 
601
 * (7)
 
602
 * t^v*bessel_j(v,a*t) -> 2^v/sqrt(%pi)*gamma(v+1/2)*a^v*r^(-2*v-1)
 
603
 *
 
604
 * Maxima doesn't recognize the relationship between gamma(2*v+1) and
 
605
 * gamma(v+1).
 
606
 */
 
607
radcan(specint(t^v*bessel_j(v,a*t)*exp(-p*t),t));
 
608
a^v*gamma(2*v+1)/((p^2+a^2)^((2*v+1)/2)*2^v*gamma(v+1))$
 
609
 
 
610
/*
 
611
 * (9)
 
612
 * t^u*bessel_j(v,a*t) -> gamma(u+v+1)*r^(-u-1)*P(u,-v,p/r)
 
613
 */
 
614
(assume(v+u+1>0),true);
 
615
true$
 
616
(assume(a>0),true);
 
617
true$
 
618
 
 
619
radcan(specint(t^u*bessel_j(v,a*t)*exp(-p*t),t));
 
620
assoc_legendre_p(-u-1,-v,p/sqrt(p^2+a^2))*gamma(v+u+1)
 
621
        /((p^2+a^2)^((u+1)/2))$
 
622
 
 
623
/*
 
624
 * (25)
 
625
 * bessel_j(0,2*sqrt(a)*sqrt(t)) -> exp(-a/p)/p
 
626
 */
 
627
specint(bessel_j(0,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 
628
%e^-(a/p)/p$
 
629
 
 
630
/*
 
631
 * (27)
 
632
 * sqrt(t)*bessel_j(1,2*sqrt(a)*sqrt(t)) -> sqrt(a)*p^(-2)*exp(-a/p)
 
633
 */
 
634
specint(sqrt(t)*bessel_j(1,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 
635
sqrt(a)*%e^-(a/p)/p^2$
 
636
 
 
637
/*
 
638
 * (29)
 
639
 * t^(-1/2)*bessel_j(1,2*sqrt(a)*sqrt(t)) ->
 
640
 *    sqrt(%pi)/sqrt(p)*exp(-a/2/p)*bessel_i(v/2,a/2/p)
 
641
 */
 
642
specint(t^(-1/2)*bessel_j(1,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 
643
sqrt(%pi)*bessel_i(1/2,a/(2*p))*%e^-(a/(2*p))/sqrt(p)$
 
644
 
 
645
/*
 
646
 * (30)
 
647
 * t^(v/2)*bessel_j(v,2*sqrt(a)*sqrt(t)) ->
 
648
 *    a^(v/2)/p^(v+1)*exp(-a/p)
 
649
 */
 
650
specint(t^(v/2)*bessel_j(v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 
651
a^(v/2)*p^(-v-1)*%e^-(a/p)$
 
652
 
 
653
/*
 
654
 * (31)
 
655
 * t^(-v/2)*bessel_j(v,2*sqrt(a)*sqrt(t)) ->
 
656
 *    exp(%i*v*%pi)*p^(v-1)/a^(v/2)/gamma(v)*exp(-a/p)*
 
657
 *     %gammagreek(v,a/p*exp(-%i*%pi)
 
658
 *
 
659
 * %gammagreek is the incomplete gamma function.
 
660
 */
 
661
specint(t^(-v/2)*bessel_j(v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 
662
p^(v-1)*%e^-(a/p)*v*%gammagreek(v,-a/p)/(a^(v/2)*(-1)^v*gamma(v+1))$
 
663
 
 
664
/*
 
665
 * (32)
 
666
 * t^(v/2-1)*bessel_j(v,2*sqrt(a)*sqrt(t)) ->
 
667
 *    a^(-v/2)*%gammagreek(v,a/p)
 
668
 */
 
669
specint(t^(v/2-1)*bessel_j(v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 
670
v*gamma(v)*%gammagreek(v,a/p)/(a^(v/2)*gamma(v+1))$
 
671
 
 
672
/*
 
673
 * (34)
 
674
 * t^(u-1/2)*bessel_j(2*v,2*sqrt(a)*sqrt(t)) ->
 
675
 *   gamma(u+v+1/2)/sqrt(a)/gamma(2*v+1)*p^(-u)*exp(-a/p/2)*%m[u,v](a/p)
 
676
 *
 
677
 * A&S 13.1.32 gives
 
678
 *
 
679
 *   %m[k,u](z) = exp(-z/2)*z^(u+1/2)*M(1/2+u-k,1+2*u,z)
 
680
 *
 
681
 * A&S 13.1.27 (Kummer Transformation):
 
682
 *
 
683
 *   M(a,b,z) = exp(z)*M(b-a,b,-z)
 
684
 *
 
685
 * So
 
686
 *
 
687
 *   %m[k,u](z) = exp(z/2)*z^(u+1/2)*M(1/2+u+k,1+2*u,-z)
 
688
 *
 
689
 * But %m[-k,u](-z) = exp(z/2)*(-z)^(u+1/2)*M(1/2+u+k,1+2*u,-z)
 
690
 *
 
691
 * Therefore
 
692
 *
 
693
 *   %m[k,u](z) = (-1)^(u+1/2)*%m[-k,u](-z)
 
694
 *
 
695
 * So the Laplace transform can also be written as
 
696
 *
 
697
 *   gamma(u+v+1/2)/sqrt(a)/gamma(2*v+1)*p^(-u)*exp(-a/p/2)
 
698
 *     *%m[-u,v](-a/p)*(-1)^(v+1/2)
 
699
 *
 
700
 * Which is the answer we produce.
 
701
 */
 
702
prefer_whittaker:true;
 
703
true$
 
704
(assume(2*v+2*u+1>0),true);
 
705
true$
 
706
 
 
707
specint(t^(u-1/2)*bessel_j(2*v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 
708
%m[-u,v](-a/p)*%e^-(a/(2*p))*(-1)^(-v-1/2)*gamma(v+u+1/2)
 
709
         /(sqrt(a)*p^u*gamma(2*v+1))$
 
710
 
 
711
/*
 
712
 * (35)
 
713
 * t^(u-1)*bessel_j(2*v,2*sqrt(a)*sqrt(t)) ->
 
714
 *    gamma(u+v)*a^v/gamma(2*v+1)/p^(u+v)*%f[1,1](u+v,2*v+1,-a/p)
 
715
 */
 
716
prefer_whittaker:false;
 
717
false$
 
718
(assume(v+u>0),true);
 
719
true$
 
720
 
 
721
specint(t^(u-1)*bessel_j(2*v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 
722
a^v*p^(-v-u)*gamma(v+u)*%f[1,1]([v+u],[2*v+1],-a/p)/gamma(2*v+1)$
 
723
 
 
724
/*
 
725
 * (45)
 
726
 * bessel_y(v,a*t) ->
 
727
 *    a^v*cot(v*%pi)/r*R^(-v)-a^(-v)*csc(v*%pi)/r*R^v
 
728
 * For |Re v| < 1.
 
729
 *
 
730
 */
 
731
expand(factor(radcan(specint(exp(-p*t)*bessel_y(1/6,a*t),t))));
 
732
sqrt(3)*a^(1/6)/(sqrt(p^2+a^2)*(sqrt(p^2+a^2)+p)^(1/6))
 
733
        -2*(sqrt(p^2+a^2)+p)^(1/6)/(a^(1/6)*sqrt(p^2+a^2))$
 
734
 
 
735
(assume(v1 > 0, v1 < 1), true);
 
736
true$
 
737
expand(factor(radcan(specint(exp(-p*t)*bessel_y(v1,a*t),t))));
 
738
a^v1*cot(%pi*v1)/(sqrt(p^2+a^2)*(sqrt(p^2+a^2)+p)^v1)
 
739
       -(sqrt(p^2+a^2)+p)^v1/(a^v1*sqrt(p^2+a^2)*sin(%pi*v1)) $
 
740
 
 
741
 
 
742
/*
 
743
 * (42)
 
744
 *
 
745
 * t^(lam-1)*bessel_j(2*u,2*sqrt(a)*sqrt(t))*bessel_j(2*v,2*sqrt(a)*sqrt(t)) ->
 
746
 *    gamma(lam+u+v)/gamma(2*u+1)/gamma(2*v+1)*a^(u+v)/p^(lam+u+v)
 
747
 *      *%f[3,3]([u+v+1/2,u+v+1,lam+u+v],[2*u+1,2*v+1,2*u+2*v+1],-4*a/p)
 
748
 *
 
749
 */
 
750
(assume(u>0,v>0,lam>0),true);
 
751
true$
 
752
specint(t^(lam-1)*bessel_j(2*u,2*sqrt(a)*sqrt(t))*bessel_j(2*v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 
753
a^(v+u)*p^(-v-u-lam)*gamma(v+u+lam)
 
754
              *%f[3,3]([v+u+1/2,v+u+1,v+u+lam],[2*u+1,2*v+1,2*v+2*u+1],-4*a/p)
 
755
        /(gamma(2*u+1)*gamma(2*v+1))$ 
 
756
 
 
757
 
 
758
/*
 
759
 * (44)
 
760
 *
 
761
 * bessel_y(0,a*t) -> -2/%pi/sqrt(p^2+a^2)*asinh(p/a)
 
762
 *
 
763
 * Maxima returns 
 
764
 *
 
765
 * -2/%pi/sqrt(p^2+a^2)*legendre_q(0,p/sqrt(p^2+a^2))
 
766
 *
 
767
 * But legendre_q(0,p/r) = log((1+p/r)/(1-p/r))/2, where r = sqrt(p^2+a^2). 
 
768
 * This simplifies to log((1+p/r)/a) = log(p/a+sqrt(1+(p/a)^2)) = asinh(p/a).
 
769
 *
 
770
 * So we have -2/%pi/sqrt(p^2+a^2)*asinh(p/a).
 
771
 */
 
772
specint(bessel_y(0,a*t)*exp(-p*t),t);
 
773
-2/%pi/sqrt(p^2+a^2)*legendre_q(0,p/sqrt(p^2+a^2)) $
 
774
 
 
775
/*
 
776
 * (46)
 
777
 *
 
778
 * t*bessel_y(0,a*t)
 
779
 *     -> 2/%pi/(p^2+a^2)*(1-p/sqrt(p^2+a^2)*log((p+sqrt(p^2+a^2))/a))
 
780
 *
 
781
 * Maxima returns
 
782
 *
 
783
 *    -2/%pi/(p^2+a^2)*legendre_q(1,p/sqrt(p^2+a^2))
 
784
 *
 
785
 * But
 
786
 *     legendre_q(1,p/r) = p/r/2*log((r+p)/(r-p)) - 1
 
787
 *                       = p/r*log((p+r)/a) - 1
 
788
 *
 
789
 * So, the transform is
 
790
 *
 
791
 *     -2/%pi/(p^2+a^2)*(p/r*log((p+r)/a) - 1)
 
792
 *
 
793
 *       = 2/%pi/(p^2+a^2)*(1-p/sqrt(p^2+a^2)*log((p+sqrt(p^2+a^2))/a))
 
794
 */
 
795
specint(t*bessel_y(0,a*t)*exp(-p*t),t);
 
796
-2/%pi/(p^2+a^2)*legendre_q(1,p/sqrt(p^2+a^2)) $
 
797
 
 
798
/*
 
799
 * (47)
 
800
 *
 
801
 * t*bessel_y(1,a*t)
 
802
 *     -> -2/%pi/(p^2+a^2)*(p/a+a/sqrt(p^2+a^2)*log((p+sqrt(p^2+a^2))/a)
 
803
 *
 
804
 * Maxima returns
 
805
 *   -4/%pi/(p^2+a^2)*assoc_legendre_q(1,-1,p/r)
 
806
 *
 
807
 * But
 
808
 *
 
809
 *   assoc_legendre_q(1,-1,z) 
 
810
 *      = sqrt(1-z^2)/2/(z^2-1)*((z^2-1)*log((1+z)/(1-z)) - 2*z)
 
811
 *
 
812
 * So
 
813
 *
 
814
 *   assoc_legendre_q(1,-1,p/r) 
 
815
 *      = (a/r)/2*(-(r/a)^2)*(-(a/r)^2*log((R/a)^2)-2*p/r)
 
816
 *      = 1/2*(p/a+a/r*log(R/a))
 
817
 *
 
818
 * where R = p + r
 
819
 *
 
820
 * Finally, the transform is
 
821
 *
 
822
 * -2/%pi/(p^2+a^2)*(p/a+a/r*log(R/a))
 
823
 *
 
824
 * as expected.
 
825
 *  
 
826
 */
 
827
specint(t*bessel_y(1,a*t)*exp(-p*t),t);
 
828
-4/%pi/(p^2+a^2)*assoc_legendre_q(1,-1,p/sqrt(p^2+a^2))$
 
829
 
 
830
/*
 
831
 * Some tests for step7
 
832
 */
 
833
 
 
834
/*
 
835
 * F(s,s+1/2;2*s+1;z) can be transformed to F(s,s+1/2;2*s+2;z) via
 
836
 * A&S 15.2.6.  And we know that F(s,s+1/2;2*s+1;z) =
 
837
 * 2^(2*s)/(1+sqrt(1-z))^(2*s).
 
838
 *
 
839
 * A&S 15.2.6 says 
 
840
 * F(a,b;c+n;z) = poch(c,n)/poch(c-a,n)/poch(c-b,n)*(1-z)^(c+n-a-b)
 
841
 *                 *diff((1-z)^(a+b-c)*F(a,b;c;z),z,n)
 
842
 *
 
843
 * F(s,s+1/2;2*s+2;z)
 
844
 *    = poch(2*s+1,1)/poch(s+1,1)/poch(s+1/2,1)*(1-z)^(3/2)
 
845
 *       *diff((1-z)^(-1/2)*F(s,s+1/2;2*s+1;z),z,1)
 
846
 *
 
847
 * 
 
848
 */
 
849
 
 
850
hgfred([s,s+1/2],[2*s+2],z)
 
851
  -(2*s+1)/(s+1)/(s+1/2)*(1-z)^(3/2)*diff((1-z)^(-1/2)
 
852
    *hgfred([s,s+1/2],[2*s+1],z),z);
 
853
0$
 
854
 
 
855
/*
 
856
 * F(s,s+1/2;2*s+1;z) can be transformed to F(s+2,s+1/2;2*s+1;z) via
 
857
 * A&S 15.2.3:
 
858
 *    F(a+n,b;c;z) = z^(1-a)/poch(a,n)*diff(z^(a+n-1)*F(a,b;c;z),z,n)
 
859
 *
 
860
 * F(s+2,s+1/2;2*s+1,z) 
 
861
 *    = z^(1-s)/s/(s+1)*diff(z^(s+1)*F(s,s+1/2;2*s+1;z),z,2)
 
862
 */
 
863
 
 
864
hgfred([s+2,s+1/2],[2*s+1],z)
 
865
  - z^(1-s)/s/(s+1)*diff(z^(s+1)*hgfred([s,s+1/2],[2*s+1],z),z,2);
 
866
0$
 
867
 
 
868
/* Tests for Airy functions */
 
869
closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
 
870
closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
 
871
 
 
872
/* Derivatives of Airy functions */
 
873
diff(airy_ai(x),x);
 
874
airy_dai(x);
 
875
diff(airy_dai(x),x);
 
876
x*airy_ai(x);
 
877
diff(airy_bi(x),x);
 
878
airy_dbi(x);
 
879
diff(airy_dbi(x),x);
 
880
x*airy_bi(x);
 
881
 
 
882
/* A&S 10.4.1  Airy functions satisfy the Airy differential equation */
 
883
diff(airy_ai(x),x,2)-x*airy_ai(x);
 
884
0;
 
885
diff(airy_bi(x),x,2)-x*airy_bi(x);
 
886
0;
 
887
 
 
888
/* A&S 10.4.4  Normalization of Airy Ai function */
 
889
(c1:3^(-2/3)/gamma(2/3), closeto(airy_ai(0)-c1,1.0e-15));
 
890
true;
 
891
closeto(airy_bi(0)/sqrt(3)-c1,1.0e-15);
 
892
true;
 
893
 
 
894
/* A&S 10.4.5  Normalization of Airy Bi function */
 
895
(c2:3^(-1/3)/gamma(1/3),closeto(-airy_dai(0)-c2,1.0e-15));
 
896
true;
 
897
closeto(airy_dbi(0)/sqrt(3)-c2,1.0e-14);
 
898
true;
 
899
 
 
900
/* A&S 10.4.10 - Wronskian */
 
901
AS_10_4_10(z):=expand(airy_ai(z)*airy_dbi(z)-airy_bi(z)*airy_dai(z)-1/%pi);
 
902
AS_10_4_10(z):=expand(airy_ai(z)*airy_dbi(z)-airy_bi(z)*airy_dai(z)-1/%pi);
 
903
closeto(AS_10_4_10(1),1.0e-15);
 
904
true;
 
905
closeto(AS_10_4_10(1+%i),1.0e-15);
 
906
true;
 
907
closeto(AS_10_4_10(%i),1.0e-15);
 
908
true;
 
909
closeto(AS_10_4_10(-1+%i),2.0e-15);
 
910
true;
 
911
closeto(AS_10_4_10(-1),1.0e-15);
 
912
true;
 
913
closeto(AS_10_4_10(-1-%i),2.0e-15);
 
914
true;
 
915
closeto(AS_10_4_10(-%i),1.0e-15);
 
916
true;
 
917
closeto(AS_10_4_10(1-%i),1.0e-15);
 
918
true;
 
919
 
 
920
/* A&S 10.4.14 - only for z>0 ? */
 
921
AS_10_4_14(z):=block([y:(2/3)*(z)^(3/2)],airy_ai(z)-(sqrt(z/3)*bessel_k(1/3,y)/%pi));
 
922
AS_10_4_14(z):=block([y:(2/3)*(z)^(3/2)],airy_ai(z)-(sqrt(z/3)*bessel_k(1/3,y)/%pi));
 
923
closeto(AS_10_4_14(1),1.0e-15);
 
924
true;
 
925
closeto(AS_10_4_14(2),1.0e-15);
 
926
true;
 
927
closeto(AS_10_4_14(5),1.0e-15);
 
928
true;
 
929
closeto(AS_10_4_14(10),1.0e-15);
 
930
true;
 
931
 
 
932
/* A&S 10.4.15 - only for z<0 ? */
 
933
AS_10_4_15(z):=block([y:(2/3)*(-z)^(3/2)],airy_ai(z)-(1/3)*sqrt(-z)*(bessel_j(1/3,y)+bessel_j(-1/3,y)));
 
934
AS_10_4_15(z):=block([y:(2/3)*(-z)^(3/2)],airy_ai(z)-(1/3)*sqrt(-z)*(bessel_j(1/3,y)+bessel_j(-1/3,y)));
 
935
closeto(AS_10_4_15(-1),1.0e-15);
 
936
true;
 
937
closeto(AS_10_4_15(-2),1.0e-15);
 
938
true;
 
939
closeto(AS_10_4_15(-5),1.0e-14);
 
940
true;
 
941
closeto(AS_10_4_15(-10),1.0e-15);
 
942
true;
 
943
 
 
944
/* A&S 10.4.16 - only for z>0 ? */
 
945
AS_10_4_16(z):=block([y:(2/3)*(z)^(3/2)],airy_dai(z)+(z/sqrt(3))*bessel_k(2/3,y)/%pi);
 
946
AS_10_4_16(z):=block([y:(2/3)*(z)^(3/2)],airy_dai(z)+(z/sqrt(3))*bessel_k(2/3,y)/%pi);
 
947
closeto(AS_10_4_16(1),1.0e-15);
 
948
true;
 
949
closeto(AS_10_4_16(2),1.0e-15);
 
950
true;
 
951
closeto(AS_10_4_16(5),1.0e-15);
 
952
true;
 
953
closeto(AS_10_4_16(10),1.0e-15);
 
954
true;
 
955
 
 
956
/* A&S 10.4.17 - only for z<0 ?, Appears to be a sign error in A&S */
 
957
AS_10_4_17(z):=block([y:(2/3)*(-z)^(3/2)],airy_dai(z)-(z/3)*(bessel_j(-2/3,y)-bessel_j(2/3,y)));
 
958
AS_10_4_17(z):=block([y:(2/3)*(-z)^(3/2)],airy_dai(z)-(z/3)*(bessel_j(-2/3,y)-bessel_j(2/3,y)));
 
959
closeto(AS_10_4_17(-1),1.0e-15);
 
960
true;
 
961
closeto(AS_10_4_17(-2),1.0e-15);
 
962
true;
 
963
closeto(AS_10_4_17(-5),1.0e-14);
 
964
true;
 
965
closeto(AS_10_4_17(-10),1.0e-15);
 
966
true;
 
967
 
 
968
/* Test that complex float arguments are evaluated */
 
969
airy_ai(%i);
 
970
airy_ai(%i);
 
971
floatnump(realpart(airy_ai(1.0*%i)));
 
972
true;
 
973
 
 
974
kill(c1,c2,AS_10_4_10,AS_10_4_14,AS_10_4_15,AS_10_4_16,AS_10_4_17);
 
975
done;
 
976
 
 
977
/* End of Airy function tests */
 
978
 
 
979
/* Numerical tests of gamma function. */
 
980
 
 
981
/* A&S Table 1.1, to 15 DP */
 
982
closeto(gamma(1/2)-1.772453850905516,2e-15);
 
983
true;
 
984
closeto(gamma(1/3)-2.678938534707748,2e-15);
 
985
true;
 
986
closeto(gamma(7/4)-0.919062526848883,1e-15);
 
987
true;
 
988
 
 
989
/* Complex values.  Checked against A&S Table 6.7 to 12 DP */
 
990
closeto(gamma(1+%i)-(0.49801566811836-0.15494982830181*%i),1e-14);
 
991
true;
 
992
closeto(gamma(1+5*%i)-(-0.00169966449436-0.00135851941753*%i),1e-14);
 
993
true;
 
994
closeto(gamma(2+3*%i)-(-0.08239527266561+0.09177428743526*%i),1e-14);
 
995
true;
 
996
 
 
997
/* Test numerical evaluation of Bessel functions
 
998
 * When the order is 0, and the arg is a float, we should produce a number.
 
999
 */
 
1000
closeto(bessel_j(0,1.0) - .7651976865579666, 1e-14);
 
1001
true;
 
1002
 
 
1003
closeto(bessel_y(0,1.0) - .08825696421567691, 1e-14);
 
1004
true;
 
1005
 
 
1006
closeto(bessel_i(0,1.0) - 1.266065877752009, 1e-14);
 
1007
true;
 
1008
 
 
1009
closeto(bessel_k(0,1.0) - .4210244382407085, 1e-14);
 
1010
true;
 
1011
 
 
1012
kill(closeto);
 
1013
done;