1jn(3,4);
2jn(3,4); /*0.1320341839226408 ;*/
3j0(1);
4j0(1); /*0.7651976865579665 ;*/
5bessel_j(0,1.0);
60.7651976865579665 ;
7
8/* BESSEL is gone. This makes sure it's gone. */
9bessel(2,3);
10bessel(2,3)$
11
12bessel(2,3.0);
13bessel(2, 3.0);
14
15bessel_j(3,2.0);
160.1289432494744021;
17
18bessel_j(-5,x);
19-bessel_j(5,x);
20bessel_j(-n,x);
21bessel_j(-n,x);
22
23diff(bessel_j(0,x),x);
24-bessel_j(1,x);
25diff(bessel_j(1,x),x);
26(bessel_j(0,x)-bessel_j(2,x))/2;
27diff(bessel_y(0,x),x);
28-bessel_y(1,x);
29diff(bessel_y(1,x),x);
30(bessel_y(0,x)-bessel_y(2,x))/2;
31diff(bessel_i(0,x),x);
32bessel_i(1,x);
33diff(bessel_i(1,x),x);
34(bessel_i(2,x)+bessel_i(0,x))/2;
35diff(bessel_k(0,x),x);
36-bessel_k(1,x);
37diff(bessel_k(1,x),x);
38-((bessel_k(2,x)+bessel_k(0,x))/2);
39
40/* Maxima used to get this wrong: A&S 9.1.65 */
41diff(bessel_y(v,z),v);
42cot(%pi*v)*('diff(bessel_j(v,z),v,1)-%pi*bessel_y(v,z))-csc(%pi*v)*'diff(bessel_j(-v,z),v,1)-%pi*bessel_j(v,z);
43
44besselexpand:true;
45true$
46
47bessel_j(1/2,x);
48sqrt(2/%pi)*sin(x)/sqrt(x);
49bessel_j(-1/2,x);
50sqrt(2/%pi)*cos(x)/sqrt(x);
51
52bessel_j(3/2,x);
53(sqrt(2)*sqrt(x)*(sin(x)/x^2-cos(x)/x))/sqrt(%pi);
54
55bessel_j(5/2,x);
56(sqrt(2)*sqrt(x)*((3/x^3-1/x)*sin(x)-(3*cos(x))/x^2))/sqrt(%pi)$
57
58ratsimp(bessel_j(-3/2,x) - (2*(-1/2)/x*bessel_j(-1/2,x)-bessel_j(1/2,x)));
590$
60
61ratsimp(bessel_j(-5/2,x) - (2*(-3/2)/x*bessel_j(-3/2,x)-bessel_j(-1/2,x)));
620$
63
64ratsimp(bessel_y(1/2,x) + sqrt(2/%pi)*cos(x)/sqrt(x));
650$
66
67ratsimp(bessel_y(-1/2,x) - sqrt(2/%pi)*sin(x)/sqrt(x));
680$
69
70ratsimp(bessel_y(3/2,x) - (2*(1/2)/x*bessel_y(1/2,x)-bessel_y(-1/2,x)));
710$
72
73ratsimp(bessel_y(5/2,x) - (2*(3/2)/x*bessel_y(3/2,x)-bessel_y(1/2,x)));
740$
75
76ratsimp(bessel_y(-3/2,x) - (2*(-1/2)/x*bessel_y(-1/2,x)-bessel_y(1/2,x)));
770$
78
79ratsimp(bessel_y(-5/2,x) - (2*(-3/2)/x*bessel_y(-3/2,x)-bessel_y(-1/2,x)));
800$
81
82bessel_i(1/2,x);
83sqrt(2*x/%pi)*sinh(x)/x;
84
85bessel_i(-1/2,x);
86sqrt(2*x/%pi)*cosh(x)/x;
87
88ratsimp(bessel_i(3/2,x) - (-2*(1/2)/x*bessel_i(1/2,x)+bessel_i(-1/2,x)));
890$
90
91ratsimp(bessel_i(5/2,x) -(-2*(3/2)/x*bessel_i(3/2,x)+bessel_i(1/2,x)));
920$
93
94ratsimp(bessel_i(-3/2,x) -(2*(-1/2)/x*bessel_i(-1/2,x)+bessel_i(1/2,x)));
950$
96
97ratsimp(bessel_i(-5/2,x) - (2*(-3/2)/x*bessel_i(-3/2,x)+bessel_i(-1/2,x)));
980$
99
100ratsimp(bessel_k(1/2,x) - sqrt(%pi/2)*%e^(-x)/sqrt(x));
1010$
102
103ratsimp(bessel_k(-1/2,x)- sqrt(%pi/2)*%e^(-x)/sqrt(x));
1040$
105
106ratsimp(bessel_k(3/2,x) - (2*(1/2)/x*bessel_k(1/2,x)+bessel_k(-1/2,x)));
1070$
108
109ratsimp(bessel_k(5/2,x) -(2*(3/2)/x*bessel_k(3/2,x)+bessel_k(1/2,x)));
1100$
111
112ratsimp(bessel_k(-3/2,x) -(-2*(-1/2)/x*bessel_k(-1/2,x)+bessel_k(1/2,x)));
1130$
114
115ratsimp(bessel_k(-5/2,x) -(-2*(-3/2)/x*bessel_k(-3/2,x)+bessel_k(-1/2,x)));
1160$
117
118(assume(p>0),true);
119true$
120(assume(4*p+a>0),true);
121true$
122besselexpand:false;
123false$
124
125specint(t^(1/2)*%e^(-a*t/4)*%e^(-p*t),t);
126sqrt(%pi)/(2*(p+a/4)^(3/2));
127
128prefer_whittaker:true;
129true$
130
131/*
132 * Reference:  Table of Integral Transforms.
133 */
134
135/*
136 * f24p146:
137 *
138 * t^(v-1)*exp(-t^2/(8*a))
139 *     -> gamma(v)*2^v*a^(v/2)*exp(a*p^2)*D[-v](2*p*sqrt(a))
140 *
141 * We have v = 7/4, a = b/4 so the result should be
142 *
143 *   gamma(7/4)*2^(7/4)*(b/4)^(7/8)*exp(b*p^2/4)*D[-7/4](p*sqrt(b))
144 *     = 3/4*gamma(3/4)*b^(7/8)*exp(b*p^2/4)*D[-7/4](p*sqrt(b))
145 *
146 * But
147 *
148 *   D[v](z) = 2^(v/2+1/4)*z^(-1/2)*%w[v/2+1/4,1/4](z^2/2)
149 *
150 * and
151 *
152 *   %w[k,u](z) = gamma(-2*u)/gamma(1/2-u-k)*%m[k,u](z)
153 *                 + gamma(2*u)/gamma(1/2+u-k)*%m[k,-u](z)
154 *
155 * Thus,
156 *
157 *   D[-7/4](p*sqrt(b)) = %w[-5/8,1/4](b*p^2/2)/(2^(5/8)*b^(1/4)*sqrt(p))
158 *
159 * and
160 *
161 *   %w[-5/8,1/4](b*p^2)
162 *     = 8/3/gamma(3/8)*%m[-5/8,-1/4](b*p^2/2)
163 *         - 2*sqrt(%pi)/gamma(7/8)*%m[-5/8,1/4](b*p^2/2)
164 *
165 * And finally the transform is
166 *
167 * 3*gamma(3/4)*b^(5/8)*exp(b*p^2/4)/(4*2^(5/8)*sqrt(p))
168 *   *(8*sqrt(%pi)/3/gamma(3/8)*%m[-5/8,-1/4](b*p^2/2)
169 *      - 2*sqrt(%pi)/gamma(7/8)*%m[-5/8,1/4](b*p^2/2))
170 */
171(assume(b>0),true);
172true$
173
174specint(t^(3/4)*%e^(-t^2/2/b)*%e^(-p*t),t);
175
176/*
177-sqrt(%pi)*b^(5/8)
178          *(3*gamma(3/8)*gamma(3/4)*%e^(b*p^2/4)*%m[-5/8,1/4](b*p^2/2)
179           -4*gamma(3/4)*gamma(7/8)*%e^(b*p^2/4)
180             *%m[-5/8,-1/4](b*p^2/2))
181 /(2*2^(5/8)*gamma(3/8)*gamma(7/8)*sqrt(p))$
182 */
183
1843*gamma(3/4)*b^(7/8)
185 *%e^(b*p^2/4)
186 *(2^(19/8)*sqrt(%pi)*%m[-5/8,-1/4](b*p^2/2)/(3*gamma(3/8)*b^(1/4)*sqrt(p))
187  -2^(3/8)*sqrt(%pi)*%m[-5/8,1/4](b*p^2/2)/(gamma(7/8)*b^(1/4)*sqrt(p)))
188 /4;
189
190/*
191 * Sec. 4.5, formula (33):
192 *
193 * t^(-1/2)*exp(-2*sqrt(a)*sqrt(t)) ->
194 *    sqrt(%pi)/sqrt(p)*exp(a/p)*erfc(sqrt(a)/sqrt(p))
195 */
196ratsimp(specint(t^(-1/2)*%e^(-2*a^(1/2)*t^(1/2))*%e^(-p*t),t));
197-sqrt(%pi)*(erf(sqrt(a)/sqrt(p))-1)*%e^(a/p)/sqrt(p)$
198
199/*
200 * The Laplace transform of sin(a*t)*cosh(b*t^2) can be derived by
201 * expressing sin and cosh in terms of exponential functions.  We end
202 * up with terms of the form:
203 *
204 *   exp(+/-%i*a*t)*exp(+/-b*t^2)
205 *
206 * All of these can be computed using formula 24, p. 146 of Tables of
207 * Integral Transforms, which handle functions of the form
208 * t^(v-1)*exp(-t^2/8/a).
209 *
210 * But, we have terms of the form exp(b*t^2-p*t+%i*a*t).  I don't
211 * think this converges, so the Laplace transform doesn't exist if b >
212 * 0.
213 *
214 */
215/*
216radcan(specint(sin(a*t)*cosh(b*t^2)*%e^(-p*t),t));
217-%e^-((p^2+2*%i*a*p+a^2)/(4*b))*(sqrt(%pi)*%e^((2*%i*a*p+a^2)/(2*b))*erf((%i
218 *p+a)/(2*sqrt(b)))-sqrt(%pi)*%e^(a^2/(2*b))*erf((%i*p-a)/(2*sqrt(b)))+sqrt
219 (%pi)*%i*%e^((p^2+2*%i*a*p)/(2*b))*erf((p+%i*a)/(2*sqrt(b)))-sqrt(%pi)*%i*%e
220 ^(p^2/(2*b))*erf((p-%i*a)/(2*sqrt(b)))+(sqrt(%pi)*%i-sqrt(%pi)*%i*%e^(%i*a
221 *p/b))*%e^(p^2/(2*b))-sqrt(%pi)*%e^((2*%i*a*p+a^2)/(2*b))+sqrt(%pi)*%e^(a^2/
222 (2*b)))/(8*sqrt(b)) $
223*/
224
225/*
226 * Sec 4.14, formula (27):
227 *
228 * t^(1/2)*bessel_j(1,2*a^(1/2)*t^(1/2)) ->
229 *    sqrt(a)/p^2*exp(-a/p)
230 */
231
232specint(t^(1/2)*bessel_j(1,2*a^(1/2)*t^(1/2))*%e^(-p*t),t);
233sqrt(a)*%e^-(a/p)/p^2$
234
235/*
236 * Sec 4.14, formula (3):
237 *
238 * t^2*bessel_j(v,a*t) ->
239 *    ((v^2-1)/r^3 + 3*p*(p+v*r)/r^5)*(a/R)^v
240 *
241 * where r = sqrt(p^2+a^2), R = p + r;
242 *
243 * (Maxima can't currently compute this transform for general v due to a bug
244 * in hyp.lisp.)
245 * This bug is no longer present after correction of legf24 in hyp.lisp.
246 */
247factor(ratsimp(specint(t^2*bessel_j(1,a*t)*%e^(-p*t),t)));
2483*a*p/(p^2+a^2)^(5/2) $
249
250(/* This is the Laplace transform of the Struve H function, see
251  http://dlmf.nist.gov/Draft/ST/about_ST.8.13.html */
252 2/(%pi*p)-2*p*log(p/(sqrt(p^2+1)-1))/(%pi*sqrt(p^2+1)),
253 /* And this should be the same as the specint of the next test below */
254 -diff(%%,p),
255 ev(fullratsimp(%%),logexpand:all));
256-((sqrt(p^2+1)*(2*p^2*log(sqrt(p^2+1)-1)-2*p^2*log(p))-2*p^2-2) /(%pi*p^6+2*%pi*p^4+%pi*p^2))$
257
258(ev(fullratsimp(specint(t*struve_h(1,t)*%e^(-p*t),t)),logexpand:all),
259 ratsimp(%%/%));
2601$
261
262/*
263 *
264 * From the comments for hstf in hypgeo.lisp:
265 *
266 * struve_h(1,t) = 2/sqrt(%pi)*(t/2)^2/gamma(1+3/2)*%f[1,2]([1],[3/2,5/2],-t^2/4)
267 *
268 * So
269 *
270 * struve_h(1,sqrt(t)) = 2/(3*%pi)*t*%f[1,2]([1],[3/2,5/2],-t/4)
271 *
272 * and the integrand is
273 *
274 * 2/(3*%pi)*t^(5/2)*exp(-p*t)*%f[1,2]([1],[3/2,5/2],-t/4).
275 *
276 * From the f19p220, the Laplace transform of this, with s = 7/2,
277 * c=-1/4, k = 1, is
278 *
279 * 2/(3*%pi)*gamma(7/2)/p^(7/2)*%f[2,2]([1,7/2],[3/2,5/2],-1/4/p)
280 *
281 * From the derivation of SPLITPFQ, we can simplify this
282 * hypergeometric function.
283 *
284 * %f[2,2]([1,7/2],[3/2,5/2],z) =
285 *
286 *      1
287 *     sum z^k/poch(5/2,k)*binomial(1,k) *diff(%f[2,2]([1,5/2],[3/2,5/2],z,k)
288 *     k=0
289 *
290 * But %f[2,2]([1,5/2],[3/2,5/2],z) = %f[1,1]([1],[3/2],z)
291 * and Maxima knows how to compute this.
292 */
293ratsimp(specint(t^(3/2)*struve_h(1,t^(1/2))*%e^(-p*t),t));
294-%e^-(1/(4*p))*(sqrt(%pi)*sqrt(p)
295				*(8*%i*erf(%i/(2*sqrt(p)))*p
296				 -%i*erf(%i/(2*sqrt(p))))
297		      -2*p*%e^(1/(4*p)))
298	/(8*sqrt(%pi)*p^(9/2)) $
299
300/* Trivial result because %ibes is not bessel_i
301   After specializing the pattern match of arbpow1 we get a more
302   correct noun form. The result is adjusted. DK.
303 */
304specint(t*%ibes[0](a*t/2)*%ibes[1](a*t/2)*%e^(-p*t),t);
305/* %ibes[0](a*t/2)*%ibes[1](a*t/2)/p^2 $ */
306specint(t*%ibes[0](a*t/2)*%ibes[1](a*t/2)*%e^(-p*t),t);
307
308/*
309 * t^(3/4)*bessel_j(1/2,t)*bessel_j(1/4,t)
310 *
311 * Luke gives
312 *
313 * bessel_j(u,t)*bessel_j(v,t)
314 *    = (z/2)^(u+v)/gamma(u+1)/gamma(v+1)
315 *        * %f[2,3]([(u+v+1)/2,(u+v+2)/2],[u+1,v+1,u+v+1],-z^2)
316 *
317 * So the integrand is
318 *
319 * 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)
320 *
321 * f19p220 gives
322 *
323 * t^(3/2)*%f[2,3]([7/8,11/8],[3/2,5/4,7/4],-t^2)
324 *
325 *    -> gamma(5/2)*p^(-5/2)
326 *         *%f[4,3]([7/8,11/8,5/2/2,(5/2+1)/2],[3/2,5/4,7/4],-4/p^2)
327 *    =  gamma(5/2)*p^(-5/2)
328 *         $%f[2,1]([7/8,11/8],[3/2],-4/p^2)
329 *
330 * And we know %f[2,1]([7/8,11/8],[3/2],z) is
331 *
332 * -2*(1/(sqrt(z)+1)^(3/4)-1/(1-sqrt(z))^(3/4))/(3*sqrt(z))
333 *
334 * Applying all of this gives the expected answer below.
335 */
336
337specint(t^(3/4)*bessel_j(1/2,t)*bessel_j(1/4,t)*%e^(-p*t),t);
338(2^(1/4)*%i*(1/((2*%i)/p+1)^(3/4)-1/(1-(2*%i)/p)^(3/4)))/(gamma(1/4)*p^(3/2))$
339
340/*
341 * Not sure this is right.  We can convert bessel_y to bessel_j,
342 * multiply them together and use the results for products of bessel_j
343 * functions.
344 *
345 * bessel_y(1/2,sqrt(t)) = -bessel_j(-1/2,sqrt(t))
346 *
347 * And maxima should be able to compute the transform of
348 *
349 * t^(5/2)*bessel_j(-1/2,sqrt(t))^2
350 *
351 * Or note that bessel_y(1/2,sqrt(t)) =
352 * -sqrt(2/%pi)*cos(sqrt(t))/t^(1/4).  Then the integrand becomes
353 *
354 * 2/%pi*t^2*cos(sqrt(t))^2
355 *
356 * And maxima should know how to compute the transform of this.
357 *
358 * Unfortunately, the transforms of these two approaches don't agree.
359 * Yuck!
360 * After revision 1.65 of hypgeo.lisp it works as expected.
361 */
362result:factor(ratsimp(specint(t^(5/2)*bessel_y(1/2,t^(1/2))^2*%e^(-p*t),t)));
363%e^-(1/p)*(16*p^3*%e^(1/p)-18*p^2*%e^(1/p)+4*p*%e^(1/p)
364                                +15*sqrt(%pi)*%i*erf(%i/sqrt(p))*p^(5/2)
365                                -20*sqrt(%pi)*%i*erf(%i/sqrt(p))*p^(3/2)
366                                +4*sqrt(%pi)*%i*erf(%i/sqrt(p))*sqrt(p))
367       /(4*%pi*p^6);
368
369/* This is equal to the Laplace transform of 2/%pi*t^2*cos(sqt(t))^2 */
370expand(result-specint(t^(5/2)*bessel_y(1/2,t^(1/2))^2*%e^(-p*t),t)),
371   besselexpand:true;
372 0;
373
374/*
375 * See formula (42), p. 187:
376 *
377 * t^(lam-1)*bessel_j(2*u,2*sqrt(a)*sqrt(t))*bessel_j(2*v,2*sqrt(a)*sqrt(t))
378 *
379 *    -> 2*gamma(lam+u+v)*a^(u+v)/gamma(2*u+1)/gamma(2*v+1)/p^(lam+u+v)
380 *        *%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)
381 *
382 * with Re(lam + u + v) > 0.
383 *
384 * So, we have lam = 3/2, u=v=1/4, a = 1/4, we get
385 *
386 * 4/%pi/p^2*%f[3,3]([1,3/2,2],[3/2,3/2,2],-1/p)
387 *  = 4/%pi/p^2*%f[1,1]([1],[3/2],-1/p)
388 *
389 * And %f[1,1]([1],[3/2],-1/p) is
390 *
391 *     -sqrt(%pi)*%i*erf(%i*sqrt(1/p))*%e^-(1/p)/(2*sqrt(1/p))
392 *
393 * So, the final result is:
394 *
395 * -2*%i*erf(%i/sqrt(p))*%e^-(1/p)/(sqrt(%pi)*p^(3/2))
396 *
397 * But we also have
398 *
399 * bessel_j(u,t)*bessel_j(v,t)
400 *    = (z/2)^(u+v)/gamma(u+1)/gamma(v+1)
401 *        * %f[2,3]([(u+v+1)/2,(u+v+2)/2],[u+1,v+1,u+v+1],-z^2)
402 *
403 * So bessel_j(1/2,sqrt(t))^2 is
404 *
405 *    2/%pi*%f[2,3]([1,3/2],[3/2,3/2,2],-t)*sqrt(t)
406 *
407 * So the integrand is
408 *
409 *    2/%pi*t*%f[2,3]([1,3/2],[3/2,3/2,2],-t)
410 *     = 2/%pi*t*%f[1,2]([1],[3/2,2],-t)
411 *
412 * f19p220 then gives us the desired transform:
413 *
414 *    t*%f[1,2]([1],[3/2,2],-t)
415 *      -> gamma(2)*p^(-2)*%f[2,2]([1,2],[3/2,2],-1/p)
416 *
417 *      = p^(-2)*%f[1,1]([1],[3/2],-1/p)
418 *
419 * So the final answer is
420 *
421 *    -%i*erf(%i/sqrt(p))*%e^-(1/p)/(sqrt(%pi)*p^(3/2))
422 *
423 * Hmm.  This differs from formula 42 above.  I think there's a bug in
424 * formula 42, and it should be divided by 2.
425 *
426 * If we use the expression for the product of Bessel functions and
427 * f19p220, we can easily derive the result of formula 42, except,
428 * we're missing the factor of 2.  So, I think formula 42 is wrong.
429 */
430
431specint(t^(1/2)*bessel_j(1/2,t^(1/2))^2*%e^(-p*t),t);
432-%i*erf(%i/sqrt(p))*%e^-(1/p)/(sqrt(%pi)*p^(3/2)) $
433
434/*
435 * See formula (8), section 4.16 of Table of Integral Transforms:
436 *
437 * t^u*bessel_i(v,a*t) -> gamma(u+v+1)*s^(-u-1)*assoc_legendre_p(u,-v,p/s)
438 *
439 * where s = sqrt(p^2-a^2);
440 */
441factor(ratsimp(specint(t^(1/2)*bessel_i(1,t)*%e^(-p*t),t)));
4423*sqrt(%pi)*assoc_legendre_p(1/2,-1,p/sqrt(p^2-1))/(4*(p^2-1)^(3/4))$
443
444/*
445 * hankel_1(2/3,sqrt(t)) = bessel_j(2/3,sqrt(t))+%i*bessel_y(2/3,sqrt(t))
446 *
447 * Formula (34) below gives:
448 *
449 * t^(u-1/2)*bessel_j(2*v,2*sqrt(a)*sqrt(t)) ->
450 *   gamma(u+v+1/2)/sqrt(a)/gamma(2*v+1)*p^(-u)*exp(-a/p/2)*%m[u,v](a/p)
451 *
452 * Formula (50) gives
453 *
454 * t^(u-1/2)*bessel_y(2*v,2*sqrt(a)*sqrt(t)) ->
455 *   1/sqrt(a)*p^(-u)*exp(-a/p/2)*
456 *     (tan((u-v)*%pi)*gamma(u+v-1/2)/gamma(2*v+1) * %m[u,v](a/p)
457 *       - sec((u-v)*%pi)*%w[u,v](a/p))
458 *
459 * But A&S 13.1.34 says
460 *
461 * %w[k,u](z) = gamma(-2*u)/gamma(1/2-u-k)*%m[k,u](z)
462 *               + gamma(2*u)/gamma(1/2+u-k)*%m[k,-u](z)
463 *
464 */
465
466ratsimp(specint(t*hankel_1(2/3,t^(1/2))*%e^(-p*t),t));
467/* Because of revision 1.110 of hyp.lisp Maxima knows in addition
468 *    hgfred([7/3],[5/3],-1/(4*x)),
469 * the result is in terms of the bessel_i function.
470
471-4*%i*gamma(1/3)*%m[-3/2,1/3](-1/(4*p))*%e^-(1/(8*p))/(3*(-1)^(5/6)*sqrt(3)
472 *gamma(2/3)*p^(3/2))+4*gamma(1/3)*%m[-3/2,1/3](-1/(4*p))*%e^-(1/(8*p))/(3*(
473 -1)^(5/6)*gamma(2/3)*p^(3/2))-8*%i*gamma(2/3)*%m[-3/2,-1/3](-1/(4*p))*%e^-
474 (1/(8*p))/(3*(-1)^(1/6)*sqrt(3)*gamma(1/3)*p^(3/2)) $ */
475
476(((-1)^(1/6)*2^(2/3)*sqrt(3)*%i-3*(-1)^(1/6)*2^(2/3))
477        *gamma(1/3)^2*gamma(5/6)*bessel_i(11/6,-1/(8*p))
478        +10*(-1)^(5/6)*sqrt(3)*4^(2/3)*%i*gamma(1/6)*gamma(2/3)^2
479           *bessel_i(7/6,-1/(8*p))
480        +(45*(-1)^(1/6)*2^(2/3)-5*(-1)^(1/6)*2^(2/3)*3^(3/2)*%i)
481         *gamma(1/3)^2*gamma(5/6)*bessel_i(5/6,-1/(8*p))
482        +((9*(-1)^(1/6)*2^(5/3)-(-1)^(1/6)*2^(5/3)*3^(3/2)*%i)
483         *bessel_i(-1/6,-1/(8*p))
484         +(5*(-1)^(1/6)*2^(5/3)*sqrt(3)*%i-15*(-1)^(1/6)*2^(5/3))
485          *bessel_i(-7/6,-1/(8*p)))
486         *gamma(1/3)^2*gamma(5/6)
487        +(((-1)^(5/6)*sqrt(3)*4^(2/3)*%i*bessel_i(-11/6,-1/(8*p))
488         -5*(-1)^(5/6)*3^(3/2)*4^(2/3)*%i*bessel_i(-5/6,-1/(8*p)))
489         *gamma(1/6)
490         -2*(-1)^(5/6)*3^(3/2)*4^(2/3)*%i*gamma(1/6)*bessel_i(1/6,-1/(8*p)))
491         *gamma(2/3)^2)
492        *%e^-(1/(8*p))
493        /(15*2^(13/3)*4^(1/3)*gamma(1/3)*gamma(2/3)*p^(7/2))/-1;
494
495/*
496 * hankel_2(3/4,t) = bessel_j(3/4,t)-%i*bessel_y(3/4,t)
497 *
498 * Sec 4.14, formula (9):
499 *
500 * t^u*bessel_j(v,a*t) -> gamma(u+v+1)*r^(-u-1)*assoc_legendre_p(u,-v,p/r)
501 *
502 * where r = sqrt(p^2+a^2)
503 *
504 * Sec 4.14, formula (48)
505 *
506 * t^u*bessel_y(v,a*t)
507 *    -> r^(-u-1)*(gamma(u+v+1)*cot(v*%pi)*assoc_legendre_p(u,-v,p/r)
508 *                  -gamma(u-v+1)*csc(v*%pi)*assoc_legendre_p(u,v,p/r))
509 *
510 * So,
511 *
512 * t^(1/2)*bessej_j(3/4,t)
513 *    -> gamma(9/4)*r^(-3/2)*assoc_legendre_p(1/2,-3/4,p/r)
514 *    =  5*gamma(1/4)/16*r^(-3/2)*assoc_legendre_p(1/2,-3/4,p/r)
515 *
516 * t^(1/2)*bessel_y(3/4,t)
517 *    -> r^(-3/2)*(gamma(9/4)*cot(3/4*%pi)*assoc_legendre_p(1/2,-3/4,p/r)
518 *                  -gamma(3/4)*csc(3/4*%pi)*assoc_legendre_p(1/2,3/4,p/r))
519 *    =  r^(-3/2)*(-gamma(9/4)*assoc_legendre_p(1/2,-3/4,p/r)
520 *                  -gamma(3/4)*sqrt(2)*assoc_legendre_p(1/2,3/4,p/r))
521 */
522
523ratsimp(specint(t^(1/2)*hankel_2(3/4,t)*%e^(-p*t),t));
524''(ratsimp(5*gamma(1/4)/16/(p^2+1)^(3/4)*assoc_legendre_p(1/2,-3/4,p/sqrt(p^2+1))
525+sqrt(2)*%i*assoc_legendre_p(1/2,3/4,p/sqrt(p^2+1))*gamma(3/4)
526       /(p^2+1)^(3/4)
527       +5*%i*gamma(1/4)*assoc_legendre_p(1/2,-3/4,p/sqrt(p^2+1))
528   /(16*(p^2+1)^(3/4))))$
529
530/*
531 * hankel_1(1/2,t) = bessel_j(1/2,t)+%i*bessel_y(1/2,t)
532 *
533 * So,
534 *
535 * t^(3/2)*bessel_j(1/2,t)
536 *     -> gamma(3/2+1/2+1)*r^(-5/2)*assoc_legendre_p(3/2,-1/2,p/r)
537 *     =  2*r^(-5/2)*assoc_legendre_p(3/2,-1/2,p/r)
538 * t^(3/2)*bessel_y(1/2,t)
539 *     -> r^(-5/2)*(gamma(3/2+1/2+1)*cot(%pi/2)*assoc_legendre_p(3/2,-1/2,p/r)
540 *                  -gamma(3/2-1/2+1)*csc(%pi/2)*assoc_legendre_p(3/2,1/2,p/r))
541 *     =  -r^(-5/2)*assoc_legendre_p(3/2,1/2,p/r))
542 *
543 * assoc_legendre_p(3/2,+/-1/2,z) can be expressed in terms of
544 * hypergeometric functions (A&S 8.1.2).
545 *
546 * assoc_legendre_p(3/2,-1/2,z)
547 *   = 1/gamma(3/2)*((z+1)/(z-1))^(-1/4)*F(-3/2,5/2;3/2;(1-z)/2)
548 *   = sqrt(2)*(z-1)^(1/4)*z*(z+1)^(1/4)/sqrt(%pi)
549 *
550 * assoc_legendre_p(3/2,1/2,z)
551 *   = 1/gamma(-1/2)*((z+1)/(z-1))^(1/4)*F(-3/2,5/2;-1/2;(1-z)/2)
552 *   = sqrt(2)*z*(2*z^2-3)/(sqrt(%pi)*(z-1)^(1/4)*(z+1)^(5/4))
553 *
554 *
555 * So the result should be
556 *
557 * t^(3/2)*bessel_j(1/2,t)
558 *    -> 4*p/(sqrt(2)*sqrt(%pi)*(p^2+1)^2)
559 *
560 * t^(3/2)*bessel_y(1/2,t)
561 *    -> -sqrt(2)*(p-1)*(p+1)/(sqrt(%pi)*(p^2+1)^2)
562 */
563
564ratsimp(specint(t^(3/2)*hankel_1(1/2,t)*%e^(-p*t),t));
565-((sqrt(%pi)*(sqrt(2)*%i*p^2-2*sqrt(2)*p-sqrt(2)*%i))/(%pi*p^4+2*%pi*p^2+%pi))$
566
567/*
568 * Formula 2, p 105:
569 *
570 * t^(u-1)*bessel_y(v,a*t)
571 *    -> -2/%pi*gamma(u+v)*(a^2+p^2)^(-u/2)
572 *         *assoc_legendre_q(u-1,-v,p/sqrt(a^2+p^2))
573 *
574 * for a > 0, Re u > |Re v|
575 *
576 * We have u = 5/2, v = 1, so the result is
577 *
578 *    -4/%pi/(p^2+a^2)*assoc_legendre_q(3/2,-1,p/sqrt(p^2+a^2))
579 *
580 * The expected result is not correct.
581 * With gamma(1+5/2) = 15*sqrt(%pi)/8 we get:
582 *
583 * 15/(4*sqrt(%pi))*(p^2+a^2)^(-5/4)*assoc_legendre_q(3/3,1,p/sqrt(p^2+a^2))
584 *
585 * That is the result of Maxima too. The example is correct.
586 */
587factor(specint(t^(5/2-1)*bessel_y(1,a*t)*%e^(-p*t),t));
588-15*assoc_legendre_q(3/2,-1,p/sqrt(p^2+a^2))/(4*sqrt(%pi)*(p^2+a^2)^(5/4));
589
590/*
591 * A&S 13.1.32:
592 *
593 *   %m[k,u](t) = exp(-t/2)*t^(u+1/2)*M(1/2+u-k,1+2*u,t)
594 *
595 * So
596 *
597 *   %m[1/2,1](t) = exp(-t/2)*t^(3/2)*M(1,3,t)
598 *
599 * and the integrand is:
600 *
601 *   t^(3/2)*%m[1/2,1](t)*exp(-p*t)
602 *      = t^3*M(1,3,t)*exp(-(p+1/2)*t)
603 *      = t^3*M(1,3,t)*exp(-p'*t)
604 *
605 * f19p220 will give us the Laplace transform of t^3*M(1,3,t):
606 *
607 *   gamma(4)/p'^4*F(1,4;3;1/p')
608 *
609 * But p' = p+1/2, so the final result is
610 *
611 *   32*(6*p-1)/(2*p-1)^2/(2*p+1)^3
612 *
613 */
614ratsimp(specint( t^(3/2)*%m[1/2,1](t)*%e^(-p*t),t));
615''(ratsimp(32*(6*p-1)/(2*p-1)^2/(2*p+1)^3));
616
617(assume(p>a),true);
618true;
619/*
620 * exp(a*t)*t^2*erf(sqrt(t))*exp(-p*t)
621 *
622 * A&S 7.1.21 gives erf(z) = 2/sqrt(%pi)*z*M(1/2,3/2,-z^2) so
623 * erf(sqrt(t)) = 2/sqrt(%pi)*sqrt(t)*M(1/2,3/2,-t)
624 *
625 * Therefore, the integrand, with p' = p-a, is
626 *
627 *   2/sqrt(%pi)*t^(5/2)*M(1/2,3/2,-t)*exp(-p'*t)
628 *
629 * Applying f19p220, the Laplace transform is
630 *
631 *   2/sqrt(%pi)*gamma(7/2)/p'^(7/2)*F(1/2,7/2;3/2;-1/p')
632 *
633 * Maxima can compute F(1/2,7/2;3/2;-1/p') and substituting p'=p-a
634 * gives us the desired answer.
635 *
636 * 15*(1/sqrt(1/(p-a)+1)-2/(3*(p-a)*(1/(p-a)+1)^(3/2))
637 *                      +1/(5*(p-a)^2*(1/(p-a)+1)^(5/2)))
638 *  /(4*(p-a)^(7/2))
639 */
640specint(%e^(a*t)*t^2*erf(t^(1/2))*%e^(-p*t),t);
64115*(1/sqrt(1/(p-a)+1)-2/(3*(p-a)*(1/(p-a)+1)^(3/2))
642                     +1/(5*(p-a)^2*(1/(p-a)+1)^(5/2)))
643 /(4*(p-a)^(7/2)) $
644
645/*
646 * Laplace transforms from Tables of Integral Transforms
647 */
648
649/*
650 * p 182, (1)
651 *
652 * bessel_j(v,a*t) -> r^(-1)*(a/R)^v
653 *
654 * where r = sqrt(p^2+a^2) and R = p + r
655 */
656(assume(v>0),true);
657true$
658
659radcan(specint(bessel_j(v,a*t)*exp(-p*t),t));
660a^v/(sqrt(p^2+a^2)*(sqrt(p^2+a^2)+p)^v)$
661
662/*
663 * (5)
664 * bessel_j(v,a*t)/t -> v^(-1)*(a/R)^v
665 *
666 * (Maxima doesn't recognize that gamma(v)/gamma(v+1) is 1/v.)
667 */
668radcan(specint(bessel_j(v,a*t)/t*exp(-p*t),t));
669a^v*gamma(v)/((sqrt(p^2+a^2)+p)^v*gamma(v+1))$
670
671/*
672 * (7)
673 * t^v*bessel_j(v,a*t) -> 2^v/sqrt(%pi)*gamma(v+1/2)*a^v*r^(-2*v-1)
674 *
675 * Maxima doesn't recognize the relationship between gamma(2*v+1) and
676 * gamma(v+1).
677 */
678radcan(specint(t^v*bessel_j(v,a*t)*exp(-p*t),t));
679a^v*gamma(2*v+1)/((p^2+a^2)^((2*v+1)/2)*2^v*gamma(v+1))$
680
681/*
682 * (9)
683 * t^u*bessel_j(v,a*t) -> gamma(u+v+1)*r^(-u-1)*P(u,-v,p/r)
684 */
685(assume(v+u+1>0),true);
686true$
687(assume(a>0),true);
688true$
689
690radcan(specint(t^u*bessel_j(v,a*t)*exp(-p*t),t));
691/* This is not the correct answer: see the formula above which is correct.
692   We had a bug in the routine legf24 in hyp.lisp.
693assoc_legendre_p(-u-1,-v,p/sqrt(p^2+a^2))*gamma(v+u+1)
694	/((p^2+a^2)^((u+1)/2))$
695 */
696
697assoc_legendre_p(u,-v,p/sqrt(p^2+a^2))*gamma(v+u+1)
698   /((p^2+a^2)^((u+1)/2))$
699
700/*
701 * (25)
702 * bessel_j(0,2*sqrt(a)*sqrt(t)) -> exp(-a/p)/p
703 */
704specint(bessel_j(0,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
705%e^-(a/p)/p$
706
707/*
708 * (27)
709 * sqrt(t)*bessel_j(1,2*sqrt(a)*sqrt(t)) -> sqrt(a)*p^(-2)*exp(-a/p)
710 */
711specint(sqrt(t)*bessel_j(1,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
712sqrt(a)*%e^-(a/p)/p^2$
713
714/*
715 * (29)
716 * t^(-1/2)*bessel_j(1,2*sqrt(a)*sqrt(t)) ->
717 *    sqrt(%pi)/sqrt(p)*exp(-a/2/p)*bessel_i(v/2,a/2/p)
718 */
719specint(t^(-1/2)*bessel_j(1,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
720sqrt(%pi)*bessel_i(1/2,a/(2*p))*%e^-(a/(2*p))/sqrt(p)$
721
722/*
723 * (30)
724 * t^(v/2)*bessel_j(v,2*sqrt(a)*sqrt(t)) ->
725 *    a^(v/2)/p^(v+1)*exp(-a/p)
726 */
727specint(t^(v/2)*bessel_j(v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
728a^(v/2)*p^(-v-1)*%e^-(a/p)$
729
730/*
731 * (31)
732 * t^(-v/2)*bessel_j(v,2*sqrt(a)*sqrt(t)) ->
733 *    exp(%i*v*%pi)*p^(v-1)/a^(v/2)/gamma(v)*exp(-a/p)*
734 *     gamma_incomplete_lower(v,a/p*exp(-%i*%pi)
735 */
736specint(t^(-v/2)*bessel_j(v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
737p^(v-1)*%e^-(a/p)*v*gamma_incomplete_lower(v,-a/p)/(a^(v/2)*(-1)^v*gamma(v+1))$
738
739/*
740 * (32)
741 * t^(v/2-1)*bessel_j(v,2*sqrt(a)*sqrt(t)) ->
742 *    a^(-v/2)*gamma_incomplete_lower(v,a/p)
743 */
744specint(t^(v/2-1)*bessel_j(v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
745v*gamma(v)*gamma_incomplete_lower(v,a/p)/(a^(v/2)*gamma(v+1))$
746
747/*
748 * (34)
749 * t^(u-1/2)*bessel_j(2*v,2*sqrt(a)*sqrt(t)) ->
750 *   gamma(u+v+1/2)/sqrt(a)/gamma(2*v+1)*p^(-u)*exp(-a/p/2)*%m[u,v](a/p)
751 *
752 * A&S 13.1.32 gives
753 *
754 *   %m[k,u](z) = exp(-z/2)*z^(u+1/2)*M(1/2+u-k,1+2*u,z)
755 *
756 * A&S 13.1.27 (Kummer Transformation):
757 *
758 *   M(a,b,z) = exp(z)*M(b-a,b,-z)
759 *
760 * So
761 *
762 *   %m[k,u](z) = exp(z/2)*z^(u+1/2)*M(1/2+u+k,1+2*u,-z)
763 *
764 * But %m[-k,u](-z) = exp(z/2)*(-z)^(u+1/2)*M(1/2+u+k,1+2*u,-z)
765 *
766 * Therefore
767 *
768 *   %m[k,u](z) = (-1)^(u+1/2)*%m[-k,u](-z)
769 *
770 * So the Laplace transform can also be written as
771 *
772 *   gamma(u+v+1/2)/sqrt(a)/gamma(2*v+1)*p^(-u)*exp(-a/p/2)
773 *     *%m[-u,v](-a/p)*(-1)^(v+1/2)
774 *
775 * Which is the answer we produce.
776 */
777prefer_whittaker:true;
778true$
779(assume(2*v+2*u+1>0),true);
780true$
781
782specint(t^(u-1/2)*bessel_j(2*v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
783%m[-u,v](-a/p)*%e^-(a/(2*p))*(-1)^(-v-1/2)*gamma(v+u+1/2)
784	 /(sqrt(a)*p^u*gamma(2*v+1))$
785
786/*
787 * (35)
788 * t^(u-1)*bessel_j(2*v,2*sqrt(a)*sqrt(t)) ->
789 *    gamma(u+v)*a^v/gamma(2*v+1)/p^(u+v)*%f[1,1](u+v,2*v+1,-a/p)
790 */
791prefer_whittaker:false;
792false$
793(assume(v+u>0),true);
794true$
795
796specint(t^(u-1)*bessel_j(2*v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
797a^v*p^(-v-u)*gamma(v+u)*%f[1,1]([v+u],[2*v+1],-a/p)/gamma(2*v+1)$
798
799/*
800 * (45)
801 * bessel_y(v,a*t) ->
802 *    a^v*cot(v*%pi)/r*R^(-v)-a^(-v)*csc(v*%pi)/r*R^v
803 * For |Re v| < 1.
804 *
805 */
806expand(factor(radcan(specint(exp(-p*t)*bessel_y(1/6,a*t),t))));
807sqrt(3)*a^(1/6)/(sqrt(p^2+a^2)*(sqrt(p^2+a^2)+p)^(1/6))
808	-2*(sqrt(p^2+a^2)+p)^(1/6)/(a^(1/6)*sqrt(p^2+a^2))$
809
810(assume(v1 > 0, v1 < 1), true);
811true$
812expand(factor(radcan(specint(exp(-p*t)*bessel_y(v1,a*t),t))));
813a^v1*cot(%pi*v1)/(sqrt(p^2+a^2)*(sqrt(p^2+a^2)+p)^v1)
814       -(sqrt(p^2+a^2)+p)^v1/(a^v1*sqrt(p^2+a^2)*sin(%pi*v1)) $
815
816
817/*
818 * (42)
819 *
820 * t^(lam-1)*bessel_j(2*u,2*sqrt(a)*sqrt(t))*bessel_j(2*v,2*sqrt(a)*sqrt(t)) ->
821 *    gamma(lam+u+v)/gamma(2*u+1)/gamma(2*v+1)*a^(u+v)/p^(lam+u+v)
822 *      *%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)
823 *
824 */
825(assume(u>0,v>0,lam>0),true);
826true$
827specint(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);
828a^(v+u)*p^(-v-u-lam)*gamma(v+u+lam)
829	      *%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)
830	/(gamma(2*u+1)*gamma(2*v+1))$
831
832
833/*
834 * (44)
835 *
836 * bessel_y(0,a*t) -> -2/%pi/sqrt(p^2+a^2)*asinh(p/a)
837 *
838 * Maxima returns
839 *
840 * -2/%pi/sqrt(p^2+a^2)*legendre_q(0,p/sqrt(p^2+a^2))
841 *
842 * But legendre_q(0,p/r) = log((1+p/r)/(1-p/r))/2, where r = sqrt(p^2+a^2).
843 * This simplifies to log((1+p/r)/a) = log(p/a+sqrt(1+(p/a)^2)) = asinh(p/a).
844 *
845 * So we have -2/%pi/sqrt(p^2+a^2)*asinh(p/a).
846 *
847 * With revision 1.64 of hypgeo.lisp we simplify the Legendre Q function.
848 * The result is equivalent to the above formula.
849 */
850specint(bessel_y(0,a*t)*exp(-p*t),t);
851
852/*-2/%pi/sqrt(p^2+a^2)*legendre_q(0,p/sqrt(p^2+a^2)) $*/
853
854 -log((sqrt(p^2+a^2)+p)/(sqrt(p^2+a^2)-p))/(%pi*sqrt(p^2+a^2));
855
856/*
857 * (46)
858 *
859 * t*bessel_y(0,a*t)
860 *     -> 2/%pi/(p^2+a^2)*(1-p/sqrt(p^2+a^2)*log((p+sqrt(p^2+a^2))/a))
861 *
862 * Maxima returns
863 *
864 *    -2/%pi/(p^2+a^2)*legendre_q(1,p/sqrt(p^2+a^2))
865 *
866 * But
867 *     legendre_q(1,p/r) = p/r/2*log((r+p)/(r-p)) - 1
868 *                       = p/r*log((p+r)/a) - 1
869 *
870 * So, the transform is
871 *
872 *     -2/%pi/(p^2+a^2)*(p/r*log((p+r)/a) - 1)
873 *
874 *       = 2/%pi/(p^2+a^2)*(1-p/sqrt(p^2+a^2)*log((p+sqrt(p^2+a^2))/a))
875 *
876 * The Legendre Q function simplifes accordingly.
877 */
878factor(specint(t*bessel_y(0,a*t)*exp(-p*t),t));
879/*-2/%pi/(p^2+a^2)*legendre_q(1,p/sqrt(p^2+a^2)) $*/
880
881(p*log((sqrt(p^2+a^2)+p)/(sqrt(p^2+a^2)-p))-2*sqrt(p^2+a^2))
882    /(-%pi*(p^2+a^2)^(3/2));
883
884/*
885 * (47)
886 *
887 * t*bessel_y(1,a*t)
888 *     -> -2/%pi/(p^2+a^2)*(p/a+a/sqrt(p^2+a^2)*log((p+sqrt(p^2+a^2))/a)
889 *
890 * Maxima returns
891 *   -4/%pi/(p^2+a^2)*assoc_legendre_q(1,-1,p/r)
892 *
893 * But
894 *
895 *   assoc_legendre_q(1,-1,z)
896 *      = sqrt(1-z^2)/2/(z^2-1)*((z^2-1)*log((1+z)/(1-z)) - 2*z)
897 *
898 * So
899 *
900 *   assoc_legendre_q(1,-1,p/r)
901 *      = (a/r)/2*(-(r/a)^2)*(-(a/r)^2*log((R/a)^2)-2*p/r)
902 *      = 1/2*(p/a+a/r*log(R/a))
903 *
904 * where R = p + r
905 *
906 * Finally, the transform is
907 *
908 * -2/%pi/(p^2+a^2)*(p/a+a/r*log(R/a))
909 *
910 * as expected.
911 *
912 */
913factor(specint(t*bessel_y(1,a*t)*exp(-p*t),t));
914/*-4/%pi/(p^2+a^2)*assoc_legendre_q(1,-1,p/sqrt(p^2+a^2))$*/
915
916(a^2*log((sqrt(p^2+a^2)+p)/(sqrt(p^2+a^2)-p))+2*p*sqrt(p^2+a^2))
917        /(%pi*a*(p^2+a^2)^(3/2));
918
919
920/*
921 * Some tests for step7
922 */
923
924/*
925 * F(s,s+1/2;2*s+1;z) can be transformed to F(s,s+1/2;2*s+2;z) via
926 * A&S 15.2.6.  And we know that F(s,s+1/2;2*s+1;z) =
927 * 2^(2*s)/(1+sqrt(1-z))^(2*s).
928 *
929 * A&S 15.2.6 says
930 * F(a,b;c+n;z) = poch(c,n)/poch(c-a,n)/poch(c-b,n)*(1-z)^(c+n-a-b)
931 *                 *diff((1-z)^(a+b-c)*F(a,b;c;z),z,n)
932 *
933 * F(s,s+1/2;2*s+2;z)
934 *    = poch(2*s+1,1)/poch(s+1,1)/poch(s+1/2,1)*(1-z)^(3/2)
935 *       *diff((1-z)^(-1/2)*F(s,s+1/2;2*s+1;z),z,1)
936 *
937 *
938 */
939
940hgfred([s,s+1/2],[2*s+2],z)
941  -(2*s+1)/(s+1)/(s+1/2)*(1-z)^(3/2)*diff((1-z)^(-1/2)
942    *hgfred([s,s+1/2],[2*s+1],z),z);
9430$
944
945/*
946 * F(s,s+1/2;2*s+1;z) can be transformed to F(s+2,s+1/2;2*s+1;z) via
947 * A&S 15.2.3:
948 *    F(a+n,b;c;z) = z^(1-a)/poch(a,n)*diff(z^(a+n-1)*F(a,b;c;z),z,n)
949 *
950 * F(s+2,s+1/2;2*s+1,z)
951 *    = z^(1-s)/s/(s+1)*diff(z^(s+1)*F(s,s+1/2;2*s+1;z),z,2)
952 */
953
954hgfred([s+2,s+1/2],[2*s+1],z)
955  - z^(1-s)/s/(s+1)*diff(z^(s+1)*hgfred([s,s+1/2],[2*s+1],z),z,2);
9560$
957
958/* Tests for Airy functions */
959closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
960closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
961
962/* Derivatives of Airy functions */
963diff(airy_ai(x),x);
964airy_dai(x);
965diff(airy_dai(x),x);
966x*airy_ai(x);
967diff(airy_bi(x),x);
968airy_dbi(x);
969diff(airy_dbi(x),x);
970x*airy_bi(x);
971
972/* Integrals of Airy functions */
973integrate(airy_ai(z),z);
974hypergeometric([1/3],[2/3,4/3],z^3/9)*z/(3^(2/3)*gamma(2/3))
975   -3^(1/6)*gamma(2/3)*hypergeometric([2/3],[4/3,5/3],z^3/9)*z^2/(4*%pi);
976integrate(airy_dai(z),z);
977airy_ai(z);
978integrate(airy_bi(z),z);
9793^(2/3)*gamma(2/3)*hypergeometric([2/3],[4/3,5/3],z^3/9)*z^2/(4*%pi)
980   +hypergeometric([1/3],[2/3,4/3],z^3/9)*z/(3^(1/6)*gamma(2/3));
981integrate(airy_dbi(z),z);
982airy_bi(z);
983
984/* A&S 10.4.1  Airy functions satisfy the Airy differential equation */
985diff(airy_ai(x),x,2)-x*airy_ai(x);
9860;
987diff(airy_bi(x),x,2)-x*airy_bi(x);
9880;
989
990/* A&S 10.4.4  Normalization of Airy Ai function */
991(c1:3^(-2/3)/gamma(2/3), closeto(airy_ai(0)-c1,1.0e-15));
992true;
993closeto(airy_bi(0)/sqrt(3)-c1,1.0e-15);
994true;
995
996/* A&S 10.4.5  Normalization of Airy Bi function */
997(c2:3^(-1/3)/gamma(1/3),closeto(-airy_dai(0)-c2,1.0e-15));
998true;
999closeto(airy_dbi(0)/sqrt(3)-c2,1.0e-14);
1000true;
1001
1002/* Exact values A&S 10.4.4 and 10.4.5 */
1003airy_ai(0);
10041/(3^(2/3)*gamma(2/3));
1005airy_dai(0);
1006-(1/(3^(1/3)*gamma(1/3)));
1007airy_bi(0);
10081/(3^(1/6)*gamma(2/3));
1009airy_dbi(0);
10103^(1/6)/gamma(1/3);
1011
1012/* A&S 10.4.10 - Wronskian */
1013AS_10_4_10(z):=expand(airy_ai(z)*airy_dbi(z)-airy_bi(z)*airy_dai(z)-1/%pi);
1014AS_10_4_10(z):=expand(airy_ai(z)*airy_dbi(z)-airy_bi(z)*airy_dai(z)-1/%pi);
1015closeto(AS_10_4_10(1),1.0e-15);
1016true;
1017closeto(AS_10_4_10(1+%i),1.0e-15);
1018true;
1019closeto(AS_10_4_10(%i),1.0e-15);
1020true;
1021closeto(AS_10_4_10(-1+%i),2.0e-15);
1022true;
1023closeto(AS_10_4_10(-1),1.0e-15);
1024true;
1025closeto(AS_10_4_10(-1-%i),2.0e-15);
1026true;
1027closeto(AS_10_4_10(-%i),1.0e-15);
1028true;
1029closeto(AS_10_4_10(1-%i),1.0e-15);
1030true;
1031
1032/* A&S 10.4.14 - only for z>0 ? */
1033AS_10_4_14(z):=block([y:(2/3)*(z)^(3/2)],airy_ai(z)-(sqrt(z/3)*bessel_k(1/3,y)/%pi));
1034AS_10_4_14(z):=block([y:(2/3)*(z)^(3/2)],airy_ai(z)-(sqrt(z/3)*bessel_k(1/3,y)/%pi));
1035closeto(AS_10_4_14(1),1.0e-15);
1036true;
1037closeto(AS_10_4_14(2),1.0e-15);
1038true;
1039closeto(AS_10_4_14(5),1.0e-15);
1040true;
1041closeto(AS_10_4_14(10),1.0e-15);
1042true;
1043
1044/* A&S 10.4.15 - only for z<0 ? */
1045AS_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)));
1046AS_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)));
1047closeto(AS_10_4_15(-1),1.0e-15);
1048true;
1049closeto(AS_10_4_15(-2),1.0e-15);
1050true;
1051closeto(AS_10_4_15(-5),1.0e-14);
1052true;
1053closeto(AS_10_4_15(-10),1.0e-15);
1054true;
1055
1056/* A&S 10.4.16 - only for z>0 ? */
1057AS_10_4_16(z):=block([y:(2/3)*(z)^(3/2)],airy_dai(z)+(z/sqrt(3))*bessel_k(2/3,y)/%pi);
1058AS_10_4_16(z):=block([y:(2/3)*(z)^(3/2)],airy_dai(z)+(z/sqrt(3))*bessel_k(2/3,y)/%pi);
1059closeto(AS_10_4_16(1),1.0e-15);
1060true;
1061closeto(AS_10_4_16(2),1.0e-15);
1062true;
1063closeto(AS_10_4_16(5),1.0e-15);
1064true;
1065closeto(AS_10_4_16(10),1.0e-15);
1066true;
1067
1068/* A&S 10.4.17 - only for z<0 ?, Appears to be a sign error in A&S */
1069AS_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)));
1070AS_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)));
1071closeto(AS_10_4_17(-1),1.0e-15);
1072true;
1073closeto(AS_10_4_17(-2),1.0e-15);
1074true;
1075closeto(AS_10_4_17(-5),1.0e-14);
1076true;
1077closeto(AS_10_4_17(-10),1.0e-15);
1078true;
1079
1080/* Test that complex float arguments are evaluated */
1081airy_ai(%i);
1082airy_ai(%i);
1083floatnump(realpart(airy_ai(1.0*%i)));
1084true;
1085
1086kill(c1,c2,AS_10_4_10,AS_10_4_14,AS_10_4_15,AS_10_4_16,AS_10_4_17);
1087done;
1088
1089/* End of Airy function tests */
1090
1091/* Numerical tests of gamma function. */
1092
1093/* A&S Table 1.1, to 15 DP */
1094closeto(gamma(1/2)-1.772453850905516,2e-15);
1095true;
1096closeto(gamma(1/3)-2.678938534707748,3.0e-15);
1097true;
1098closeto(gamma(7/4)-0.919062526848883,1e-15);
1099true;
1100
1101/* Complex values.  Checked against A&S Table 6.7 to 12 DP */
1102closeto(gamma(1+%i)-(0.49801566811836-0.15494982830181*%i),1e-14);
1103true;
1104closeto(gamma(1+5*%i)-(-0.00169966449436-0.00135851941753*%i),1e-14);
1105true;
1106closeto(gamma(2+3*%i)-(-0.08239527266561+0.09177428743526*%i),1e-14);
1107true;
1108
1109/* Test numerical evaluation of Bessel functions
1110 * When the order is 0, and the arg is a float, we should produce a number.
1111 */
1112closeto(bessel_j(0,1.0) - .7651976865579666, 1e-14);
1113true;
1114
1115closeto(bessel_y(0,1.0) - .08825696421567691, 1e-14);
1116true;
1117
1118closeto(bessel_i(0,1.0) - 1.266065877752009, 1e-14);
1119true;
1120
1121closeto(bessel_k(0,1.0) - .4210244382407085, 1e-14);
1122true;
1123
1124/*
1125 * Tests for failed cases to see if we're returning the noun forms
1126 */
1127/* fail-on-f24p146test */
1128specint(t^(-3/2)*exp(-t^2/8/a)*exp(-p*t),t);
1129'specint(t^(-3/2)*exp(-t^2/8/a)*exp(-p*t),t);
1130
1131/* fail-on-f35p147test
1132   Because we modified the construction of a noun form, we get a sligthly
1133   different noun form as result. DK. */
1134specint((2*t)^(-3/2)*exp(-2*sqrt(a)*sqrt(t))*exp(-p*t),t);
1135/*'specint(exp(-p*t -2*sqrt(a)*sqrt(t))/(8*t^(3/2)),t);*/
1136'specint(%e^(-p*t-2*sqrt(a)*sqrt(t))/(2*sqrt(2)*t^(3/2)),t);
1137
1138/* fail-on-f29p146test */
1139specint(t^(v-1)*exp(1/8/t)*exp(-p*t),t);
1140'specint(t^(v-1)*exp(1/8/t)*exp(-p*t),t);
1141
1142/* fail-in-arbpow (aka f1p137test) */
1143specint(t^(-1)*exp(-p*t),t);
1144'specint(exp(-p*t)/t,t);
1145
1146/* fail-in-f2p105vcond */
1147(assume(p>3),0);
11480;
1149
1150specint(8*t^1*exp(3*t)*bessel_y(3,t)*exp(-p*t),t);
1151'specint(8*bessel_y(3,t)*t*%e^(3*t-p*t),t)$
1152
1153/* fail-in-f50cond */
1154specint(8*t^1*exp(3*t)*bessel_y(8,7*sqrt(t))*exp(-p*t),t);
1155'specint(8*bessel_y(8,7*sqrt(t))*t*%e^(3*t-p*t),t);
1156
1157/* fail-in-dionarghyp-y */
1158specint(8*t^1*exp(3*t)*bessel_y(8,7*t^(3/2))*exp(-p*t),t);
1159'specint(8*bessel_y(8,7*t^(3/2))*t*%e^(3*t-p*t),t);
1160
1161/* The additionally phase factor in the calculation vanish, because of
1162   the modificaton of the transformation to Bessel J in the code. DK.
1163*/
1164(assume(t>0,v2>1), radcan(specint(bessel_i(-v2,t)*exp(-p*t),t)));
1165/*'specint((%e^((%i*%pi*v2-2*p*t)/2)*bessel_i(-v2,t))/(-1)^(v2/2),t);*/
1166'(specint(bessel_i(-v2,t)*exp(-p*t),t));
1167
1168/* Verify fix for [ 1877522 ] erf(1.0),nouns wrong; causes plot2d(erf) to fail
1169 */
1170erf(1.0), nouns;
11710.8427007929497148;
1172
1173erf(1.0), erf;
11740.8427007929497148;
1175
1176erf(1.0);
11770.8427007929497148;
1178
1179/* [ 789059 ] simpexpt problem: sign called on imag arg */
1180(-(-1)^(1/6))^(1/2);
1181sqrt(-(-1)^(1/6));
1182
1183/* Further tests of bessel functions */
1184/* The following numerical values are evaluated with the evaluation tool
1185   of the website www.functions.wolfram.com with a precision of 16 digits
1186   1998-2014 Wolfram Research, Inc. */
1187
1188(test_bessel(actual, ref, digits) := closeto(realpart(actual)-realpart(ref), 10^(-digits)) and closeto(imagpart(actual)-imagpart(ref), 10^(-digits)), 0);
11890;
1190
1191/* Numerical values for the bessel function J with negative order */
1192
1193test_bessel(bessel_j(-1,-2.0),0.5767248077568734, 15);
1194true;
1195
1196test_bessel(bessel_j(-1,2.0), -0.5767248077568734, 15);
1197true;
1198
1199test_bessel(bessel_j(-1,-1.5), 0.5579365079100996, 15);
1200true;
1201
1202test_bessel(bessel_j(-1,1.5), -0.5579365079100996, 15);
1203true;
1204
1205test_bessel(bessel_j(-1.5, -2.0), -0.3956232813587035*%i, 15);
1206true;
1207
1208test_bessel(bessel_j(-1.5, 2.0), -0.3956232813587035, 15);
1209true;
1210
1211test_bessel
1212  (bessel_j(-1.8, -1.5),- 0.2033279902093184 - 0.1477264320209275 * %i,15);
1213true;
1214
1215test_bessel(bessel_j(-1.8, 1.5), -0.251327217627129314,15);
1216true;
1217
1218test_bessel(bessel_j(-2,-1.5), 0.2320876721442147, 15);
1219true;
1220
1221test_bessel(bessel_j(-2,1.5), 0.2320876721442147, 15);
1222true;
1223
1224test_bessel(bessel_j(-2.5,-1.5), -1.315037204805194 * %i,15);
1225true;
1226
1227test_bessel(bessel_j(-2.5,1.5), 1.315037204805194,15);
1228true;
1229
1230test_bessel
1231  (bessel_j(-2.3,-1.5), 0.5949438455752484 - 0.8188699527453657 * %i,14);
1232true;
1233
1234test_bessel(bessel_j(-2.3,1.5), 1.012178926325313, 14);
1235true;
1236
1237/* Numerical values for the bessel function J with positive order */
1238
1239test_bessel(bessel_j(1.5,1.0), 0.2402978391234270, 15);
1240true;
1241
1242test_bessel(bessel_j(1.5,-1.0), -0.2402978391234270 * %i, 15);
1243true;
1244
1245test_bessel(bessel_j(1.8,1.0), 0.1564953153109239, 14);
1246true;
1247
1248test_bessel
1249  (bessel_j(1.8,-1.0), 0.1266073696266034 - 0.0919856383926216 * %i, 15);
1250true;
1251
1252test_bessel(bessel_j(2.0,1.0), 0.1149034849319005,15);
1253true;
1254
1255test_bessel(bessel_j(2.0,-1.0),0.1149034849319005,15);
1256true;
1257
1258test_bessel(bessel_j(2.5,1.0),0.04949681022847794,15);
1259true;
1260
1261test_bessel(bessel_j(2.5,-1.0),0.04949681022847794 * %i,15);
1262true;
1263
1264/* Numerical values for the bessel function J with complex arg
1265   and positive or negative order*/
1266
1267test_bessel
1268  (bessel_j(0,1.0+%i),0.9376084768060293 - 0.4965299476091221 * %i,15);
1269true;
1270
1271test_bessel
1272  (bessel_j(1,1.0+%i),0.6141603349229036 + 0.3650280288270878 * %i,15);
1273true;
1274
1275test_bessel
1276  (bessel_j(-1,1.0+%i),-0.6141603349229036 - 0.3650280288270878 * %i,14);
1277true;
1278
1279test_bessel
1280  (bessel_j(2,1.0+%i),0.0415798869439621 + 0.2473976415133063 * %i,15);
1281true;
1282
1283test_bessel
1284  (bessel_j(-2,1.0+%i),0.0415798869439621 + 0.2473976415133063 * %i,15);
1285true;
1286
1287test_bessel
1288  (bessel_j(2.3,1.0+%i),-0.0141615213034667 + 0.1677798241687935 * %i,15);
1289true;
1290
1291test_bessel
1292  (bessel_j(-2.3,1.0+%i),0.1920598664138632 - 0.5158676904105332 * %i,14);
1293true;
1294
1295/* Numerical values for the bessel function J with complex order */
1296
1297test_bessel
1298  (bessel_j(%i,1.0),1.641024179495082 - 0.437075010213683*%i,15);
1299true;
1300
1301test_bessel
1302  (bessel_j(%i,1.5),1.401883276281807 + 0.473362399311655*%i,15);
1303true;
1304
1305test_bessel
1306  (bessel_j(1.0+%i,-1.0),-0.01142279482478010 + 0.02390070064911069*%i,15);
1307true;
1308
1309test_bessel
1310  (bessel_j(1.5*%i,-2.0),0.01925195427338360 + 0.01442616961986814*%i,15);
1311true;
1312
1313/******************************************************************/
1314/* Numerical values for the bessel function Y with negative order */
1315
1316test_bessel
1317  (bessel_y(-1,-2.0),-0.1070324315409375 + 1.1534496155137468 * %i,14);
1318true;
1319
1320test_bessel(bessel_y(-1,2.0),0.1070324315409375,15);
1321true;
1322
1323test_bessel
1324  (bessel_y(-1,-1.5),-0.4123086269739113 + 1.1158730158201993 * %i,15);
1325true;
1326
1327test_bessel(bessel_y(-1,1.5),0.4123086269739113,15);
1328true;
1329
1330test_bessel(bessel_y(-1.5, -2.0),0.4912937786871623 * %i,15);
1331true;
1332
1333test_bessel(bessel_y(-1.5, 2.0),-0.4912937786871623,15);
1334true;
1335
1336test_bessel
1337  (bessel_y(-1.8, -1.5),-0.6777414340388011 + 0.0857519944018923 * %i,14);
1338true;
1339
1340test_bessel(bessel_y(-1.8, 1.5),-0.8377344836401481,14);
1341true;
1342
1343test_bessel
1344  (bessel_y(-2,-1.5),-0.9321937597629739 + 0.4641753442884295 * %i,15);
1345true;
1346
1347test_bessel(bessel_y(-2,1.5),-0.9321937597629739,14);
1348true;
1349
1350test_bessel(bessel_y(-2.5,-1.5),0.1244463597983876 * %i,14);
1351true;
1352
1353test_bessel(bessel_y(-2.5,1.5),0.1244463597983876,15);
1354true;
1355
1356test_bessel
1357  (bessel_y(-2.3,-1.5),-0.3148570865836879 + 0.7565240896444820 * %i,14);
1358true;
1359
1360test_bessel(bessel_y(-2.3,1.5),-0.5356668704355646,14);
1361true;
1362
1363/* Numerical values for the bessel function Y with positive order */
1364
1365test_bessel(bessel_y(1.5,1.0),-1.102495575160179,14);
1366true;
1367
1368test_bessel(bessel_y(1.5,-1.0),-1.102495575160179 * %i,14);
1369true;
1370
1371test_bessel(bessel_y(1.8,1.0),-1.382351995367631,14);
1372true;
1373
1374test_bessel
1375  (bessel_y(1.8,-1.0),-1.1183462564605324 - 0.5593113771009602 * %i,14);
1376true;
1377
1378test_bessel(bessel_y(2.0,1.0),-1.650682606816254,14);
1379true;
1380
1381test_bessel
1382  (bessel_y(2.0,-1.0),-1.650682606816254 + 0.229806969863801 * %i,14);
1383true;
1384
1385test_bessel(bessel_y(2.5,1.0),-2.876387857462161,14);
1386true;
1387
1388test_bessel(bessel_y(2.5,-1.0),2.876387857462161 * %i,14);
1389true;
1390
1391
1392/* Numerical values for the bessel function Y with complex arg
1393   and positive or negative order*/
1394
1395test_bessel
1396  (bessel_y(0,1.0+%i),0.4454744889360325 + 0.7101585820037345 * %i,15);
1397true;
1398
1399test_bessel
1400  (bessel_y(1,1.0+%i),-0.6576945355913452 + 0.6298010039928844 * %i,15);
1401true;
1402
1403test_bessel
1404(bessel_y(-1,1.0+%i),0.6576945355913452 - 0.6298010039928844 * %i,15);
1405true;
1406
1407test_bessel
1408  (bessel_y(2,1.0+%i),-0.4733680205344934 + 0.5773369575804951 * %i,14);
1409true;
1410
1411test_bessel
1412  (bessel_y(-2,1.0+%i),-0.4733680205344934 + 0.5773369575804951 * %i,14);
1413true;
1414
1415test_bessel
1416  (bessel_y(2.3,1.0+%i),-0.2476879981252862 + 0.7595467103431256 * %i,15);
1417true;
1418
1419test_bessel
1420  (bessel_y(-2.3,1.0+%i),-0.1570442638685963 + 0.5821870838327466 * %i,14);
1421true;
1422
1423/* Numerical values for the bessel function Y with complex order */
1424
1425test_bessel
1426  (bessel_y(%i,1.0),-0.476556612479964 - 1.505069159110387*%i,14);
1427true;
1428
1429test_bessel
1430  (bessel_y(%i,1.5),0.5161218926267688 - 1.2857405211747503*%i,14);
1431true;
1432
1433test_bessel
1434  (bessel_y(1.0+%i,-1.0),7.708594946281541 + 1.233384674244926*%i,14);
1435true;
1436
1437test_bessel
1438  (bessel_y(1.5*%i,-2.0),3.226466016458932 + 4.267260420563194*%i,13);
1439true;
1440
1441/******************************************************************/
1442/* Numerical values for the bessel function I with negative order */
1443
1444test_bessel(bessel_i(-1,-2.0),-1.590636854637329,15);
1445true;
1446
1447test_bessel(bessel_i(-1,2.0),1.590636854637329,15);
1448true;
1449
1450test_bessel(bessel_i(-1,-1.5),-0.9816664285779076,15);
1451true;
1452
1453test_bessel(bessel_i(-1,1.5),0.9816664285779076,15);
1454true;
1455
1456test_bessel(bessel_i(-1.5, -2.0),0.9849410530002364 * %i,14);
1457true;
1458
1459test_bessel(bessel_i(-1.5, 2.0),0.9849410530002364,14);
1460true;
1461
1462test_bessel
1463  (bessel_i(-1.8, -1.5),0.2026903980307014 + 0.1472631941876387 * %i,15);
1464true;
1465
1466test_bessel(bessel_i(-1.8, 1.5),0.2505391103524365,15);
1467true;
1468
1469test_bessel(bessel_i(-2,-1.5),0.3378346183356807,15);
1470true;
1471
1472test_bessel(bessel_i(-2,1.5),0.3378346183356807,15);
1473true;
1474
1475test_bessel(bessel_i(-2.5,-1.5),-0.8015666610717216*%i,14);
1476true;
1477
1478test_bessel(bessel_i(-2.5,1.5),0.8015666610717216,14);
1479true;
1480
1481test_bessel
1482  (bessel_i(-2.3,-1.5),0.3733945265830230 - 0.5139334755917659 * %i,15);
1483true;
1484
1485test_bessel(bessel_i(-2.3,1.5),0.6352567117441516,15);
1486true;
1487
1488/* Numerical values for the bessel function I with positive order */
1489
1490test_bessel(bessel_i(1.5,1.0),0.2935253263474798,15);
1491true;
1492
1493test_bessel(bessel_i(1.5,-1.0),-0.2935253263474798*%i,13);
1494true;
1495
1496test_bessel(bessel_i(1.8,1.0),0.1871011888310777,15);
1497true;
1498
1499test_bessel
1500  (bessel_i(1.8,-1.0),0.1513680414320980 - 0.1099753194812967 * %i,14);
1501true;
1502
1503test_bessel(bessel_i(2.0,1.0),0.1357476697670383,15);
1504true;
1505
1506test_bessel(bessel_i(2.0,-1.0),0.1357476697670383,15);
1507true;
1508
1509test_bessel(bessel_i(2.5,1.0),0.05709890920304825,15);
1510true;
1511
1512test_bessel(bessel_i(2.5,-1.0),0.05709890920304825 * %i,15);
1513true;
1514
1515
1516/* Numerical values for the bessel function I with complex arg
1517   and positive or negative order*/
1518
1519test_bessel
1520  (bessel_i(0,1.0+%i),0.9376084768060293 + 0.4965299476091221 * %i,15);
1521true;
1522
1523test_bessel
1524  (bessel_i(1,1.0+%i),0.3650280288270878 + 0.6141603349229036 * %i,15);
1525true;
1526
1527test_bessel
1528  (bessel_i(-1,1.0+%i),0.3650280288270878 + 0.6141603349229036 * %i,15);
1529true;
1530
1531test_bessel
1532  (bessel_i(2,1.0+%i),-0.0415798869439621 + 0.2473976415133063 * %i,15);
1533true;
1534
1535test_bessel
1536  (bessel_i(-2,1.0+%i),-0.0415798869439621 + 0.2473976415133063 * %i,15);
1537true;
1538
1539test_bessel
1540  (bessel_i(2.3,1.0+%i),-0.0635524383467825 + 0.1559221140952053 * %i,15);
1541true;
1542
1543test_bessel
1544  (bessel_i(-2.3,1.0+%i),-0.4053256245784623 - 0.3724481230406298 * %i,14);
1545true;
1546
1547/* Numerical values for the bessel function I with complex order */
1548
1549test_bessel
1550  (bessel_i(%i,1.0),1.900799675819425 - 1.063960013554441*%i,14);
1551true;
1552
1553test_bessel
1554  (bessel_i(%i,1.5),2.495473638417463 - 0.601347501920535*%i,14);
1555true;
1556
1557test_bessel
1558  (bessel_i(1.0+%i,-1.0),-0.01096313515009349 + 0.03043920114776303*%i,15);
1559true;
1560
1561test_bessel
1562  (bessel_i(1.5*%i,-2.0),0.04238259669782487 - 0.01125055344512197*%i,15);
1563true;
1564
1565/******************************************************************/
1566/* Numerical values for the bessel function K with negative order */
1567
1568test_bessel
1569  (bessel_k(-1,-2.0),-0.139865881816522-4.997133057057809*%i,14);
1570true;
1571
1572test_bessel(bessel_k(-1,2.0),0.139865881816522,14);
1573true;
1574
1575test_bessel
1576  (bessel_k(-1,-1.5),-0.277387800456844 - 3.083996040296084*%i,14);
1577true;
1578
1579test_bessel(bessel_k(-1,1.5),0.2773878004568438,15);
1580true;
1581
1582test_bessel(bessel_k(-1.5, -2.0),-3.2741902342766302*%i,13);
1583true;
1584
1585test_bessel(bessel_k(-1.5, 2.0),0.1799066579520922,15);
1586true;
1587
1588test_bessel
1589  (bessel_k(-1.8, -1.5),0.3929372194683435 - 1.0725774293000646*%i,15);
1590true;
1591
1592test_bessel(bessel_k(-1.8, 1.5),0.4856971141526263,15);
1593true;
1594
1595test_bessel
1596  (bessel_k(-2,-1.5),0.5836559632566508 - 1.0613387550916862*%i,14);
1597true;
1598
1599test_bessel(bessel_k(-2,1.5),0.5836559632566508,15);
1600true;
1601
1602test_bessel
1603  (bessel_k(-2.3,-1.5),0.465598659425186 - 1.354876241730594*%i,14);
1604true;
1605
1606test_bessel(bessel_k(-2.3,1.5),0.7921237520153218,15);
1607true;
1608
1609/* Numerical values for the bessel function K with positive order */
1610
1611test_bessel(bessel_k(1.5,1.0),0.9221370088957891,14);
1612true;
1613
1614/* kein Wert von function-site !? Ist das eine Nullstelle??? */
1615test_bessel(bessel_k(1.5,-1.0),0,13);
1616true;
1617
1618test_bessel(bessel_k(1.8,1.0),1.275527037541854,15);
1619true;
1620
1621test_bessel
1622  (bessel_k(1.8,-1.0),1.0319230501560911 + 0.1619402612577788*%i,14);
1623true;
1624
1625test_bessel(bessel_k(2.0,1.0),1.624838898635177,14);
1626true;
1627
1628test_bessel(bessel_k(2.0,-1.0),1.624838898635177 - 0.426463882082061*%i,14);
1629true;
1630
1631test_bessel(bessel_k(2.5,1.0),3.227479531135262,14);
1632true;
1633
1634test_bessel(bessel_k(2.5,-1.0),-3.406861044815549*%i,14);
1635true;
1636
1637/* Numerical values for the bessel function K with complex arg
1638   and positive or negative order*/
1639
1640test_bessel
1641  (bessel_k(0,1.0+%i),0.0801977269465178 - 0.3572774592853303*%i,15);
1642true;
1643
1644test_bessel
1645  (bessel_k(1,1.0+%i),0.0245683055237403 - 0.4597194738011894 * %i,15);
1646true;
1647
1648test_bessel
1649  (bessel_k(-1,1.0+%i),0.0245683055237403 - 0.4597194738011894 * %i,15);
1650true;
1651
1652test_bessel
1653  (bessel_k(2,1.0+%i),-0.3549534413309312 - 0.8415652386102600 * %i,13);
1654true;
1655
1656test_bessel
1657  (bessel_k(-2,1.0+%i),-0.3549534413309312 - 0.8415652386102600 * %i,13);
1658true;
1659
1660test_bessel
1661  (bessel_k(2.3,1.0+%i),-0.6635905911278042 - 1.0258894849569300 * %i,15);
1662true;
1663
1664test_bessel
1665  (bessel_k(-2.3,1.0+%i),-0.6635905911278042 - 1.0258894849569300 * %i,13);
1666true;
1667
1668/* Numerical values for the bessel function K with complex order */
1669
1670test_bessel
1671  (bessel_k(%i,1.0),0.2894280370259921,14);
1672true;
1673
1674test_bessel
1675  (bessel_k(%i,1.5),0.1635839926633096,14);
1676true;
1677
1678test_bessel
1679  (bessel_k(1.0+%i,-1.0),-9.744252766030894 - 7.494570149760043*%i,14);
1680true;
1681
1682test_bessel
1683  (bessel_k(1.5*%i,-2.0),3.93512366711118 - 14.82183468316553*%i,13);
1684true;
1685
1686kill(closeto, test_bessel);
1687done;
1688
1689/* Numerical tests of the Bessel functions using the Wronskians
1690
1691   The Wronskians combines different types of Bessel functions and
1692   Bessel functions with negative and positve order.
1693   The results are very simple. Therefore the numerical calculation of the
1694   Wronskians is a good test for the different parts of the algorithmen.
1695
1696   Based on code by Dieter Kaiser.
1697*/
1698
1699/* Test the Wronskian.  wf is the Wronskian function.  w_true is simplified Wronskian.
1700 * eps is the max absolute error allowed, and the Wronskian is tested
1701 * for values of the arg between -zlimit and zlimit.
1702 */
1703(test_wronskian(wf, w_true, eps, zlimit) :=
1704block([badpoints : [],
1705       ratprint : false,
1706       abserr : 0,
1707       maxerr : -1],
1708  for order:-1 thru 1 step 1/10 do
1709  (
1710    for z: -zlimit thru zlimit step 1 do
1711    (
1712      if notequal(z,0.0) then
1713      (
1714        result : float(rectform(wf(float(order),z))),
1715        answer : float(rectform(w_true(float(order),z))),
1716        abserr : abs(result - answer),
1717	maxerr : max(maxerr, abserr),
1718        if abserr > eps then
1719        (
1720          badpoints : cons([[order, z], result, answer, abserr], badpoints)
1721        )
1722      )
1723    )
1724  ),
1725  /*
1726   * For debugging, if there are any bad points, return the maximum error
1727   * found as the first element.
1728   */
1729  if badpoints # [] then
1730    cons(maxerr, badpoints)
1731  else
1732    badpoints
1733), 0);
17340;
1735
1736/*******************************************************************************
1737
1738   Wronskian w_jj
1739
1740   A&S 9.1.15 : J[n+1](z)*J[-n](z)+J[n](z)*J[-(n+1)](z) = -2*sin(n*pi)/(pi*z)
1741
1742*******************************************************************************/
1743
1744w_jj(n,z) := bessel_j(n+1,z)*bessel_j(-n,z) + bessel_j(n,z)*bessel_j(-n-1,z);
1745w_jj(n,z) := bessel_j(n+1,z) *bessel_j(-n,z) + bessel_j(n,z)*bessel_j(-n-1,z);
1746
1747/* Calculation of w_jj for real argument */
1748
1749test_wronskian('w_jj, lambda([n,z], -2.0*sin(n*%pi)/(z*%pi)), 1e-14, 10);
1750[];
1751
1752
1753/* Calculation of w_jj for complex argument */
1754
1755test_wronskian(lambda([n,z], expand(w_jj(n,%i*z))),
1756               lambda([n,z],-2.0*sin(n*%pi)/(%i*z*%pi)),
1757               1e-8, 10);
1758[];
1759
1760/*******************************************************************************
1761
1762   Wronskian w_jy
1763
1764   A&S 9.1.16: J[n+1](z)*Y[n](z)-J[n](z)*Y[n+1,z] = 2/(pi*z)
1765
1766*******************************************************************************/
1767
1768w_jy(n,z) := bessel_j(n+1,z)*bessel_y(n,z) - bessel_j(n,z)*bessel_y(n+1,z);
1769w_jy(n,z) := bessel_j(n+1,z)*bessel_y(n,z) - bessel_j(n,z)*bessel_y(n+1,z);
1770
1771/* Calculation of w_yj for real argument */
1772
1773test_wronskian(w_jy, lambda([n,z], 2.0/(z*%pi)), 1e-14, 10);
1774[];
1775
1776
1777/* Calculation of w_jy for complex argument */
1778
1779test_wronskian(lambda([n,z], w_jy(n,z*%i)), lambda([n,z],2.0/(z*%i*%pi)), 1e-8, 10);
1780[];
1781
1782/*******************************************************************************
1783
1784   Wronskian w_ii
1785
1786   A&S 9.6.14: I[n](z)*I[-(n+1)](z)-I[n+1](z)*I[-n](z) = -2*sin(n*pi)/(pi*z)
1787
1788*******************************************************************************/
1789
1790w_ii(n,z) := bessel_i(n,z)*bessel_i(-n-1,z) - bessel_i(n+1,z)*bessel_i(-n,z);
1791w_ii(n,z) := bessel_i(n,z)*bessel_i(-n-1,z) - bessel_i(n+1,z)*bessel_i(-n,z);
1792
1793/* Calculation of w_ii for real argument */
1794
1795test_wronskian(w_ii, lambda([n,z],-2.0*sin(n*%pi)/(z*%pi)), 1e-10, 5);
1796[];
1797
1798/* Calculation of w_ii for complex argument */
1799
1800test_wronskian(lambda([n,z], w_ii(n,z*%i)),
1801               lambda([n,z], -2.0*sin(n*%pi)/(z*%i*%pi)),
1802               1e-10, 5);
1803[];
1804
1805/*******************************************************************************
1806
1807   Test Wronskian w_ik
1808
1809   A&S 9.6.15: I[n](z)*K[n+1](z)+I[n+1](z)*K[n](z) = 1/z
1810
1811*******************************************************************************/
1812
1813w_ik(n,z) := bessel_i(n,z)*bessel_k(n+1,z) + bessel_i(n+1,z)*bessel_k(n,z);
1814w_ik(n,z) := bessel_i(n,z)*bessel_k(n+1,z) + bessel_i(n+1,z)*bessel_k(n,z);
1815
1816/* Calculation of w_ik for real argument */
1817
1818test_wronskian(w_ik, lambda([n,z], 1/z), 1e-10, 5);
1819[];
1820
1821/* Calculation of w_ik for complex argument */
1822
1823test_wronskian(lambda([n,z], w_ik(n,z*%i)), lambda([n,z], 1/(z*%i)), 1e-12, 5);
1824[];
1825
1826/*******************************************************************************
1827
1828   Test Wronskian w_h1h2
1829
1830   A&S 9.1.17: H1[v+1](z)*H2[v](z)-H1[v](z)*H2[v+1](z) = -4*%i/(%pi*z)
1831
1832*******************************************************************************/
1833
1834w_h1h2(v,z) := hankel_1(v+1,z)*hankel_2(v,z) - hankel_1(v,z)*hankel_2(v+1,z);
1835w_h1h2(v,z) := hankel_1(v+1,z)*hankel_2(v,z) - hankel_1(v,z)*hankel_2(v+1,z);
1836
1837/* Calculation of w_h1h2 for real argument */
1838
1839test_wronskian(w_h1h2, lambda([v,z], -4*%i/(%pi*z)), 1e-14, 5);
1840[];
1841
1842/* Calculation of w_h1h2 for complex argument */
1843
1844test_wronskian(lambda([v,z], w_h1h2(v,z*%i)),
1845               lambda([v,z], -4/(%pi*z)),
1846               1e-13, 5);
1847[];
1848
1849/* Calculation of w_h1h2 for complex order and argument */
1850
1851test_wronskian(lambda([v,z], w_h1h2(v*%i,z*%i)),
1852               lambda([v,z], -4/(%pi*z)),
1853               1e-10, 5);
1854[];
1855
1856/******************************************************************************
1857
1858    Integrals of Bessel functions
1859
1860*******************************************************************************/
1861
1862integrate(bessel_j(0,x),x);
1863x*(bessel_j(0,x)*(2-%pi*struve_h(1,x))+%pi*bessel_j(1,x)*struve_h(0,x))/2;
1864
1865integrate(bessel_j(1,x),x);
1866-bessel_j(0,x);
1867
1868integrate(bessel_j(2,x),x);
1869hypergeometric([3/2],[5/2,3],-x^2/4)*x^3/24;
1870
1871integrate(bessel_j(1/2,x),x);
1872 2^(3/2)*hypergeometric([3/4],[3/2,7/4],-x^2/4)*x^(3/2)/(3*sqrt(%pi));
1873
1874/* http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/21/01/01/ */
1875integrate(bessel_y(0,x),x);
1876%pi*x*(bessel_y(1,x)*struve_h(0,x)+bessel_y(0,x)*struve_h(-1,x))/2;
1877
1878integrate(bessel_y(1,x),x);
1879-bessel_y(0,x);
1880
1881integrate(bessel_y(2,x),x);
1882%pi*x*(bessel_y(1,x)*struve_h(0,x)+bessel_y(0,x)*struve_h(-1,x))/2-2*bessel_y(1,x);
1883
1884integrate(bessel_y(3,x),x);
1885-2*bessel_y(2,x)-bessel_y(0,x);
1886
1887integrate(bessel_y(10,x),x);
1888%pi*x*(bessel_y(1,x)*struve_h(0,x)+bessel_y(0,x)*struve_h(-1,x))/2
1889  -2*'sum(bessel_y(2*i+1,x),i,0,4);
1890
1891integrate(bessel_y(11,x),x);
1892-2*'sum(bessel_y(2*i,x),i,1,5)-bessel_y(0,x);
1893
1894integrate(bessel_i(0,x),x);
1895x*(bessel_i(0,x)*(%pi*struve_l(1,x)+2)-%pi*bessel_i(1,x)*struve_l(0,x))/2;
1896
1897integrate(bessel_i(1,x),x);
1898bessel_i(0,x);
1899
1900integrate(bessel_i(2,x),x);
1901hypergeometric([3/2],[5/2,3],x^2/4)*x^3/24;
1902
1903integrate(bessel_i(1/2,x),x);
19042^(3/2)*hypergeometric([3/4],[3/2,7/4],x^2/4)*x^(3/2)/(3*sqrt(%pi));
1905
1906/* http://functions.wolfram.com/Bessel-TypeFunctions/BesselK/21/01/01/ */
1907integrate(bessel_k(0,x),x);
1908%pi*x*(bessel_k(1,x)*struve_l(0,x)+bessel_k(0,x)*struve_l(-1,x))/2;
1909
1910integrate(bessel_k(1,x),x);
1911-bessel_k(0,x);
1912
1913integrate(bessel_k(7,x),x);
19142*(-bessel_k(6,x)+bessel_k(4,x)-bessel_k(2,x))+bessel_k(0,x);
1915
1916integrate(bessel_k(8,x),x);
1917%pi*(struve_l(0,x)*bessel_k(1,x)+struve_l(-1,x)*bessel_k(0,x))*x/2
1918 +2*(-bessel_k(7,x)+bessel_k(5,x)-bessel_k(3,x)+bessel_k(1,x))$
1919
1920/******************************************************************************
1921
1922    Check the handling of realpart and impagpart for special case
1923    of the order and arg of the Bessel functions.
1924
1925*******************************************************************************/
1926
1927/* Check for the Bessel J function */
1928
1929(f(n,x):=[realpart(bessel_j(n,x)),imagpart(bessel_j(n,x))],done);
1930done;
1931
1932(declare(n,integer),assume(x>0),done);
1933done;
1934
1935f(n,x);
1936[bessel_j(n,x),0];
1937f(n,-x);
1938[bessel_j(n,-x),0];
1939
1940f(n+1/2,x);
1941[bessel_j(n+1/2,x),0];
1942f(n+1/2,-x);
1943[0,-%i*bessel_j(n+1/2,-x)];
1944
1945(declare(n_even,even),declare(n_odd,odd),%iargs:false,done);
1946done;
1947
1948f(n_even,%i);
1949[bessel_j(n_even,%i),0];
1950f(n_odd,%i);
1951[0,-%i*bessel_j(n_odd,%i)];
1952
1953f(n_even,x*%i);
1954[bessel_j(n_even,x*%i),0];
1955f(n_odd,x*%i);
1956[0,-%i*bessel_j(n_odd,x*%i)];
1957
1958f(n_even,(x+1)^2*%i);
1959[bessel_j(n_even,(x+1)^2*%i),0];
1960f(n_odd,(x+1)^2*%i);
1961[0,-%i*bessel_j(n_odd,(x+1)^2*%i)];
1962
1963(declare(j,imaginary),done);
1964done;
1965
1966f(n_even,(x+1)^2*j);
1967[bessel_j(n_even,(x+1)^2*j),0];
1968f(n_odd,(x+1)^2*j);
1969[0,-%i*bessel_j(n_odd,(x+1)^2*j)];
1970
1971/* Check the handling of realpart and imagpart for the Bessel I function */
1972
1973(f(n,x):=[realpart(bessel_i(n,x)),imagpart(bessel_i(n,x))],done);
1974done;
1975
1976f(n,x);
1977[bessel_i(n,x),0];
1978f(n,-x);
1979[bessel_i(n,-x),0];
1980
1981f(n+1/2,x);
1982[bessel_i(n+1/2,x),0];
1983f(n+1/2,-x);
1984[0,-%i*bessel_i(n+1/2,-x)];
1985
1986f(n_even,%i);
1987[bessel_i(n_even,%i),0];
1988f(n_odd,%i);
1989[0,-%i*bessel_i(n_odd,%i)];
1990
1991f(n_even,x*%i);
1992[bessel_i(n_even,x*%i),0];
1993f(n_odd,x*%i);
1994[0,-%i*bessel_i(n_odd,x*%i)];
1995
1996f(n_even,(x+1)^2*%i);
1997[bessel_i(n_even,(x+1)^2*%i),0];
1998f(n_odd,(x+1)^2*%i);
1999[0,-%i*bessel_i(n_odd,(x+1)^2*%i)];
2000
2001f(n_even,(x+1)^2*j);
2002[bessel_i(n_even,(x+1)^2*j),0];
2003f(n_odd,(x+1)^2*j);
2004[0,-%i*bessel_i(n_odd,(x+1)^2*j)];
2005
2006/* Check the handling of realpart and imagpart for the Bessel K function */
2007
2008(f(n,x):=[realpart(bessel_k(n,x)),imagpart(bessel_k(n,x))],done);
2009done;
2010
2011f(n,x);
2012[bessel_k(n,x),0];
2013f(n,-x);
2014[realpart(bessel_k(n,-x)),imagpart(bessel_k(n,-x))];
2015
2016f(n+1/2,x);
2017[bessel_k(n+1/2,x),0];
2018f(n+1/2,-x);
2019[0,-%i*bessel_k(n+1/2,-x)];
2020
2021f(n_even,%i);
2022[realpart(bessel_k(n_even,%i)),imagpart(bessel_k(n_even,%i))];
2023f(n_odd,%i);
2024[realpart(bessel_k(n_odd,%i)),imagpart(bessel_k(n_odd,%i))];
2025
2026
2027/* Check the handling of realpart and imagpart for the Bessel Y function */
2028
2029(f(n,x):=[realpart(bessel_y(n,x)),imagpart(bessel_y(n,x))],done);
2030done;
2031
2032f(n,x);
2033[bessel_y(n,x),0];
2034f(n,-x);
2035[realpart(bessel_y(n,-x)),imagpart(bessel_y(n,-x))];
2036
2037f(n+1/2,x);
2038[bessel_y(n+1/2,x),0];
2039f(n+1/2,-x);
2040[0,-%i*bessel_y(n+1/2,-x)];
2041
2042f(n_even,%i);
2043[realpart(bessel_y(n_even,%i)),imagpart(bessel_y(n_even,%i))];
2044f(n_odd,%i);
2045[realpart(bessel_y(n_odd,%i)),imagpart(bessel_y(n_odd,%i))];
2046
2047/* %iargs is false, so transformation disabled */
2048bessel_j(v, %i*x);
2049bessel_j(v, %i*x);
2050
2051bessel_i(v, %i*x);
2052bessel_i(v, %i*x);
2053
2054/* Set %iargs to true to enable transformation */
2055(%iargs:true, bessel_j(v, %i*x));
2056(%i*x)^v/x^v*bessel_i(v,x);
2057
2058bessel_i(v, %i*x);
2059(%i*x)^v/(x^v)*bessel_j(v, x);
2060
2061imagpart(bessel_j(2,%i*3.0));
20620;
2063realpart(bessel_j(3,%i*3.0));
20640;
2065
2066imagpart(bessel_i(2,%i*3.0));
20670;
2068realpart(bessel_i(3,%i*3.0));
20690;
2070
2071/*******************************************************/
2072/* Check the handling of conjugate on Bessel functions */
2073/*******************************************************/
2074
2075declare([w,z],complex,n,integer);
2076done;
2077assume(nonneg>=0,notequal(nonzero,0));
2078[nonneg>=0,notequal(nonzero,0)];
2079
2080conjugate(bessel_j(w,z));
2081'(conjugate(bessel_j(w,z)));
2082
2083conjugate(bessel_y(w,z));
2084'(conjugate(bessel_y(w,z)));
2085
2086conjugate(bessel_i(w,z));
2087'(conjugate(bessel_i(w,z)));
2088
2089conjugate(bessel_k(w,z));
2090'(conjugate(bessel_k(w,z)));
2091
2092conjugate(hankel_1(w,z));
2093'(conjugate(hankel_1(w,z)));
2094
2095conjugate(hankel_2(w,z));
2096'(conjugate(hankel_2(w,z)));
2097
2098/* Bessel functions with arguments off of the negative real axis
2099 * commute with conjugate
2100 */
2101conjugate(bessel_j(z,x+%i*nonzero));
2102'(bessel_j(conjugate(z),x-%i*nonzero));
2103
2104conjugate(bessel_j(z,nonneg));
2105'(bessel_j(conjugate(z),nonneg));
2106
2107conjugate(bessel_y(z,x+%i*nonzero));
2108'(bessel_y(conjugate(z),x-%i*nonzero));
2109
2110conjugate(bessel_y(z,nonneg));
2111'(bessel_y(conjugate(z),nonneg));
2112
2113conjugate(bessel_i(z,x+%i*nonzero));
2114'(bessel_i(conjugate(z),x-%i*nonzero));
2115
2116conjugate(bessel_i(z,nonneg));
2117'(bessel_i(conjugate(z),nonneg));
2118
2119conjugate(bessel_k(z,x+%i*nonzero));
2120'(bessel_k(conjugate(z),x-%i*nonzero));
2121
2122conjugate(bessel_k(z,nonneg));
2123'(bessel_k(conjugate(z),nonneg));
2124
2125conjugate(hankel_1(z,x+%i*nonzero));
2126'(hankel_2(conjugate(z),x-%i*nonzero));
2127
2128conjugate(hankel_1(z,nonneg));
2129'(hankel_2(conjugate(z),nonneg));
2130
2131conjugate(hankel_2(z,x+%i*nonzero));
2132'(hankel_1(conjugate(z),x-%i*nonzero));
2133
2134conjugate(hankel_2(z,nonneg));
2135'(hankel_1(conjugate(z),nonneg));
2136
2137/* Bessel functions J and I of integral order commute with conjugate */
2138conjugate(bessel_j(n,z));
2139'(bessel_j(n,conjugate(z)));
2140
2141conjugate(bessel_i(n,z));
2142'(bessel_i(n,conjugate(z)));
2143
2144remove([w,z],complex,n,integer);
2145done;
2146forget(nonneg>=0,notequal(nonzero,0));
2147[nonneg>=0,notequal(nonzero,0)];
2148
2149/* mailing list 2016-01-06: "coerce-float-fun and hashed arrays" */
2150
2151(kill(a, s, t), a:make_array(hashed), a[s]:5, 0);
21520;
2153
2154a[s];
21555;
2156
2157a[s], nouns;
21585;
2159
2160a[t];
2161false;
2162
2163a[t], nouns;
2164false;
2165
2166