1 /* rtest_laplace.mac -- test cases for Laplace transform
2 * Copyright 2012 by Robert Dodier.
3 * I release this work under terms of the GNU General Public License.
7 mycontext : newcontext (),
9 /* attempt to forestall asksign queries */
10 assume (a > 0, b > 0, c > 0),
11 assume (n > 0, q > 0),
12 assume_pos_pred : lambda ([e], true),
14 assuming (e1, e2) ::= buildq ([e1, e2], block ([foo, c1 : newcontext ()], assume (e1), foo : e2, killcontext (c1), foo)),
18 /* cases handled by LAPLACE in src/laplac.lisp */
38 laplace (exp (- u), u, v);
41 laplace (exp (- 7*u), u, v);
44 laplace (exp (- a*u), u, v);
47 laplace (exp (- u - b), u, v);
50 laplace (exp (- a*u - b), u, v);
55 laplace (sin (u), u, v);
58 laplace (sin (7*u), u, v);
61 laplace (sin (a*u), u, v);
66 laplace (cos (u), u, v);
69 laplace (cos (7*u), u, v);
72 laplace (cos (a*u), u, v);
77 laplace (sinh (u), u, v);
80 laplace (sinh (7*u), u, v);
83 laplace (sinh (a*u), u, v);
88 laplace (cosh (u), u, v);
91 laplace (cosh (7*u), u, v);
94 laplace (cosh (a*u), u, v);
99 (laplace (log (u), u, v), if equal (%%, (-1/v)*(log (v) + %gamma)) then true else %%);
102 (laplace (log (7*u), u, v), ev (if equal (%%, (1/abs(7))*(-1/(v/7)*(log (v/7) + %gamma))) then true else %%, logexpand = super));
105 (laplace (log (a*u), u, v), ev (if equal (%%, (1/abs(a))*(-1/(v/a)*(log (v/a) + %gamma))) then true else %%, logexpand = super));
110 laplace ('diff (foo (u), u), u, v);
111 v * 'laplace (foo (u), u, v) - foo (0);
113 laplace ('diff (foo (u), u, 2), u, v);
114 'laplace (foo (u), u, v) * v^2 - foo (0) * v - 'at ('diff (foo (u), u, 1), u = 0);
116 laplace ('diff (foo (u), u, 3), u, v);
117 'laplace (foo (u), u, v) * v^3 - foo (0) * v^2 - 'at ('diff (foo (u), u), u = 0) * v - 'at ('diff (foo (u), u, 2), u = 0);
119 laplace ('diff (foo (u), u, n), u, v);
120 'laplace (foo (u), u, v) * v^n - 'sum ('at ('diff (foo (u), u, k), u = 0) * v^((n - 1) - k), k, 0, n - 1);
124 laplace ('integrate (foo (w), w, 0, u), u, v);
125 (1/v)*'laplace (foo (u), u, v);
129 laplace ('sum (foo (u, i), i, 0, n), u, v);
130 'sum ('laplace (foo (u, i), u, v), i, 0, n);
134 laplace (erf (u), u, v);
135 exp (v^2/4) * (1 - erf (v/2)) / v;
137 laplace (erf (7*u), u, v);
138 exp (v^2/196) * (1 - erf (v/14)) / v;
140 laplace (erf (a*u), u, v);
141 exp (v^2/(4*a^2)) * (1 - erf (v/(2*a))) / v;
143 /* inverse transform */
145 laplace ('ilt (foo (v), v, u), u, v);
150 laplace (delta (u), u, v);
153 laplace (delta (u - 7), u, v);
156 laplace (delta (u - a), u, v);
159 laplace (delta (7*u), u, v);
162 laplace (delta (a*u), u, v);
165 laplace (delta (a*u - b), u, v);
166 exp (- b*v/a) / abs (a);
168 /* cases which specint handles
169 * I haven't verified the following rules -- these tests only verify that
170 * laplace returns the result advertised in specint's comments.
173 /* Algorithm 1: Laplace transform of c*t^v*exp(-s*t+e*f) */
175 /* Algorithm 1.1: Laplace transform of c*t^v*exp(-a*t^2)
176 * t^(v-1)*exp(-t^2/8/a)
177 * -> gamma(v)*2^v*a^(v/2)*exp(a*p^2)*D[-v](2*p*sqrt(a))
180 /* Algorithm 1.2: Laplace transform of c*t^v*exp(-a*sqrt(t))
181 * (2*t)^(v-1)*exp(-2*sqrt(a)*sqrt(t))
182 * -> gamma(2*v)*p^(-v)*exp(a/p/2)*D[-2*v](sqrt(2*a/p))
185 /* Algorithm 1.3: Laplace transform of t^v*exp(1/t)
186 * t^(v-1)*exp(-a/t/4)
187 * -> 2*(a/p/4)^(v/2)*bessel_k(v, sqrt(a)*sqrt(p))
190 /* Algorithm 1.4: Laplace transform of exp(exp(-t))
192 * -> a^(-p)*gamma(p,a)
195 laplace (exp (-a * exp (-u)), u, v);
196 gamma_incomplete_lower(v,a)/a^v$
198 /* Algorithm 1.5: Laplace transform of exp(exp(t))
200 * -> a^(-p)*gamma_incomplete(-p,a)
201 * (specint code doesn't match formula; code appears to be correct)
204 laplace (exp (-a * exp (u)), u, v);
205 a^v * gamma_incomplete (-v, a);
207 /* Algorithm 2: Laplace transform of u(t)*%e^(-p*t).
208 * Laplace transform of asin(w)
209 * Laplace transform of atan(w)
210 * Laplace transform of two Bessel J functions
211 * Laplace transform of two hankel_1 functions
212 * Laplace transform of two hankel_2 functions
213 * Laplace transform of hankel_1 * hankel_2
214 * Laplace transform of two Bessel Y functions
215 * Laplace transform of two Bessel K functions
216 * Laplace transform of Bessel K and Bessel Y functions
217 * Laplace transform of Bessel I and Bessel J functions
218 * Laplace transform of Bessel I and Hankel 1 functions
219 * Laplace transform of Bessel I and Hankel 2 functions
220 * Laplace transform of Bessel Y and Bessel J functions
221 * Laplace transform of Bessel K and Bessel J functions
222 * Laplace transform of Hankel 1 and Bessel J functions
223 * Laplace transform of Hankel 2 and Bessel J functions
224 * Laplace transform of Bessel Y and Hankel 1 functions
225 * Laplace transform of Bessel Y and Hankel 2 functions
226 * Laplace transform of Bessel K and Hankel 1 functions
227 * Laplace transform of Bessel K and Hankel 2 functions
228 * Laplace transform of Bessel I and Bessel Y functions
229 * Laplace transform of Bessel I and Bessel K functions
230 * Laplace transform of Struve H function
231 * Laplace transform of Struve L function
232 * Laplace transform of little Lommel s function
233 * Laplace transform of Lommel S function
234 * Laplace transform of Bessel Y function
235 * Laplace transform of Bessel K function
236 * Laplace transform of Parabolic Cylinder function
237 * Laplace transform of Incomplete Gamma function
238 * Laplace transform of Batemann function
239 * Laplace transform of Bessel J function
240 * Laplace transform of lower incomplete Gamma function
241 * Laplace transform of Hankel 1 function
242 * Laplace transform of Hankel 2 function
243 * Laplace transform of Whittaker M function
244 * Laplace transform of Whittaker M function
245 * Laplace transform of the Generalized Laguerre function, %l[v1,v2](w)
246 * Laplace transform for the Generalized Laguerre function
247 * Laplace transform for the Laguerre function
248 * Laplace transform of Gegenbauer function
249 * Laplace transform of Chebyshev function of the first kind
250 * Laplace transform of Chebyshev function of the second kind
251 * Laplace transform for the Hermite function, hermite(index1,arg1)
252 * Laplace transform of %p[v1,v2](w), Associated Legendre P function
253 * Laplace transform of Associated Legendre P function
254 * Laplace transform of %p[v1,v2,v3](w), Jacobi function
255 * Laplace transform of Jacobi P function
256 * Laplace transform of Associated Legendre function of the second kind
257 * Laplace transform of Associated Legendre function of the second kind
258 * Laplace transform of %p[v1](w), Legendre P function
259 * Laplace transform of Legendre P function
260 * Laplace transform of Whittaker W function
261 * Laplace transform of Whittaker W function
262 * Laplace transform of square of Bessel J function
263 * Laplace transform of square of Hankel 1 function
264 * Laplace transform of square of Hankel 2 function
265 * Laplace transform of square of Bessel Y function
266 * Laplace transform of square of Bessel K function
267 * Laplace transform of two Bessel I functions
268 * Laplace transform of Bessel I. We use I[v](w)=%i^n*J[n](%i*w).
269 * Laplace transform of square of Bessel I function
270 * Laplace transform of Erf function
271 * Laplace transform of the logarithmic function.
272 * Laplace transform of Erfc function
273 * Laplace transform of expintegral_ei.
274 * Laplace transform of expintegral_e1
275 * Laplace transform of expintegral_e
276 * Laplace transform of expintegral_si
277 * Laplace transform of expintegral_shi
278 * Laplace transform of expintegral_ci
279 * Laplace transform of expintegral_chi
280 * Laplace transform of Complete elliptic integral of the first kind
281 * Laplace transform of Complete elliptic integral of the first kind
282 * Laplace transform of Complete elliptic integral of the second kind
283 * Laplace transform of Complete elliptic integral of the second kind
284 * Laplace transform of %f[v1,v2](w1,w2,w3), Hypergeometric function
285 * Laplace transform of Hypergeometric function
286 * Laplace transform of c * t^v * (a+t)^w
287 * Laplace transform of c * t^v
290 /* Algorithm 2.1: Laplace transform of c*t^v*%e(-p*t)
292 * Table of Integral Transforms
297 * -> gamma(u+1)*p^(-u-1)
300 /* Algorithm 2.2: Laplace transform of c*t^v*(1+t)^w
303 /* Algorithm 2.3: Laplace transform of the Logarithmic function
306 * -> c*gamma(v)*s^(-v)*(psi[0](v)-log(s/a))
308 * This is the formula for an expression with log(t) scaled like 1/a*F(s/a).
310 * For the following cases we have to add further algorithm:
311 * log(1+a*x), log(x+a), log(x)^2.
314 /* Algorithm 2.4: Laplace transform of the Whittaker function
316 * Test for Whittaker W function. Simplify this if possible, or
317 * convert to Whittaker M function.
319 * We have r * %w[i1,i2](a)
323 * t^(v-1)*%w[k,u](a*t)
324 * -> gamma(u+v+1/2)*gamma(v-u+1/2)*a^(u+1/2)/
325 * (gamma(v-k+1)*(p+a/2)^(u+v+1/2)
326 * *2f1(u+v+1/2,u-k+1/2;v-k+1;(p-a/2)/(p+a/2))
329 /* Algorithm 2.5: Laplace transform of bessel_k(0,a*t)
331 * The general algorithm handles the Bessel K function for an order |v|<1.
332 * but does not include the special case v=0. Return the Laplace transform:
334 * bessel_k(0,a*t) --> acosh(s/a)/sqrt(s^2-a^2)
337 laplace (bessel_k (0, a*u), u, v);
338 acosh (v / a) / sqrt (v^2 - a^2);
340 /* Algorithm 3: Laplace transform of a hypergeometric function
341 * Table of Laplace transforms, p 220, formula 19:
343 * If m + k <= n + 1, and Re(s) > 0, the Laplace transform of
345 * t^(s-1)*%f[m,n]([a1,...,am],[p1,...,pn],(c*t)^k)
348 * gamma(s)/p^s*%f[m+k,n]([a1,...,am,s/k,(s+1)/k,...,(s+k-1)/k],[p1,...,pm],(k*c/p)^k)
350 * with Re(p) > 0 if m + k <= n, Re(p+k*c*exp(2*%pi*%i*r/k)) > 0 for r
351 * = 0, 1,...,k-1, if m + k = n + 1.
354 /* Algorithm 4: SPECIAL HANDLING OF Bessel Y for an integer order
356 * This is called for one Bessel Y function, when the order is an integer.
359 /* Algorithm 4.1: Laplace transform of t^n*bessel_y(v,a*t)
360 * v is an integer and n>=v
362 * Table of Integral Transforms
364 * Volume 2, p 105, formula 2 is a formula for the Y-transform of
366 * f(x) = x^(u-3/2)*exp(-a*x)
368 * where the Y-transform is defined by
370 * integrate(f(x)*bessel_y(v,x*y)*sqrt(x*y), x, 0, inf)
374 * -2/%pi*gamma(u+v)*sqrt(y)*(y^2+a^2)^(-u/2)
375 * *assoc_legendre_q(u-1,-v,a/sqrt(y^2+a^2))
377 * with a > 0, Re u > |Re v|.
379 * In particular, with a slight change of notation, we have
381 * integrate(x^(u-1)*exp(-p*x)*bessel_y(v,a*x)*sqrt(a), x, 0, inf)
383 * which is the Laplace transform of x^(u-1/2)*bessel_y(v,x).
385 * Thus, the Laplace transform is
387 * -2/%pi*gamma(u+v)*sqrt(a)*(a^2+p^2)^(-u/2)
388 * *assoc_legendre_q(u-1,-v,p/sqrt(a^2+p^2))
391 /* Algorithm 4.2: Laplace transform of t^n*bessel_y(v, a*sqrt(t))
393 * Table of Integral Transforms
395 * p. 188, formula 50:
397 * t^(u-1/2)*bessel_y(2*v,2*sqrt(a)*sqrt(t))
398 * -> a^(-1/2)*p^(-u)*exp(-a/2/p)
399 * * [tan((u-v)*%pi)*gamma(u+v+1/2)/gamma(2*v+1)*M[u,v](a/p)
400 * -sec((u-v)*%pi)*W[u,v](a/p)]
403 /* other examples which might or might not be special cases of stuff above
404 * collected from various tables:
407 /* http://en.wikipedia.org/wiki/Laplace_transform */
411 laplace (a * foo (u) + b * bar (u) + c * baz (u), u, v);
412 a * 'laplace (foo (u), u, v) + b * 'laplace (bar (u), u, v) + c * 'laplace (baz (u), u, v);
414 /* frequency differentiation */
416 laplace (u * foo (u), u, v);
417 - 'diff ('laplace (foo (u), u, v), v);
419 laplace (u^2 * foo (u), u, v);
420 'diff ('laplace (foo (u), u, v), v, 2);
422 laplace (u^3 * foo (u), u, v);
423 - 'diff ('laplace (foo (u), u, v), v, 3);
425 laplace (u^n * foo (u), u, v);
426 (- 1)^n * 'diff ('laplace (foo (u), u, v), v, n);
428 /* frequency integration */
430 laplace (foo (u) / u, u, v);
431 'integrate ('laplace (foo (u), u, v), u, v, inf);
435 laplace (foo (a * u), u, v);
436 (1 / a) * 'at ('laplace (foo (u), u, v), v = v / a);
438 /* frequency shifting */
440 assuming (v - a > 0, laplace (exp (a * u) * foo (u), u, v));
441 'laplace (foo (u), u, v - a);
443 (laplace (exp (a * u) * cosh (u), u, v), if equal (%%, (v - a) / ((v - a)^2 - 1)) then true else %%);
448 laplace (foo (u - a) * unit_step (u - a), u, v);
449 exp (- a * v) * 'laplace (foo (u), u, v);
451 laplace (sinh (u - a) * unit_step (u - a), u, v);
452 exp (-a * v) / (v^2 - 1);
456 laplace ('integrate (foo (w) * bar (u - w), w, 0, u), u, v);
457 'laplace (bar (u), u, v) * 'laplace (foo (u), u, v);
459 (laplace ('integrate (erf (w) * sinh (u - w), w, 0, u), u, v), if equal (%%, laplace (erf (u), u, v) * laplace (sinh (u), u, v)) then true else %%);
464 declare ([u, v], complex);
467 laplace (conjugate (foo (u)), u, v);
468 conjugate ('at ('laplace (foo (u), u, v), v = conjugate (v)));
470 /* erf is entire => conjugate commutes with erf */
472 laplace (conjugate (erf (u)), u, v);
473 exp (v^2/4) * (1 - erf (v/2)) / v;
475 /* sqrt and log are not entire => conjugate does not commute */
477 laplace (conjugate (sqrt (u)), u, v);
478 sqrt (%pi) * conjugate (1 / conjugate (v) ^ (3/2)) / 2;
480 laplace (conjugate (log (u)), u, v);
481 (- conjugate (log (conjugate (v))) - %gamma) / v;
483 remove ([u, v], complex);
488 laplace (unit_step (u), u, v);
491 laplace (unit_step (u - a), u, v);
496 laplace (u * unit_step (u), u, v);
501 makefact (laplace (u^n, u, v));
504 makefact (laplace ((u - a)^n * unit_step (u - a), u, v));
505 n! * exp (- a*v) / v^(n + 1);
507 makefact (laplace (u^n * exp (- b*u), u, v));
508 n! / (v + b)^(n + 1);
510 makefact (ratsimp (laplace ((u - a)^n * exp (- b*(u - a)) * unit_step (u - a), u, v)));
511 n!/(v + b)^(n + 1) * exp (- a*v);
516 gamma (q + 1) / v^(q + 1);
520 laplace (u^(1/n), u, v);
521 gamma (1/n + 1) / v^(1/n + 1);
523 (laplace (u^(1/n) * exp(-b*u), u, v),
524 if ratsimp (%% - gamma (1/n + 1) / (v + b)^(1/n + 1)) = 0 then true else %%);
527 /* exponential approach */
529 (laplace (1 - exp (- a*u), u, v), if equal (%%, a/(v * (v + a))) then true else %%);
534 (laplace (exp (- a*u) * sin (b*u), u, v), if equal (%%, b / ((v + a)^2 + b^2)) then true else %%);
537 (laplace (exp (- a*u) * cos (b*u), u, v), if equal (%%, (v + a) / ((v + a)^2 + b^2)) then true else %%);
540 /* Bessel function */
543 foo : laplace (bessel_j (n, a*u), u, v),
544 bar : (sqrt (v^2 + a^2) - v)^n / (a^n * sqrt (v^2 + a^2)),
545 makelist (makelist (ev (ratsimp (foo - bar), a = aa, n = nn), nn, 0, 4), aa, [1/5, 1/3, 1/2, 1, 2, 3, 5]),
546 map (lambda ([e], zeroequiv (e, v)), flatten (%%)),
550 /* http://tutorial.math.lamar.edu/Classes/DE/Laplace_Table.aspx */
552 /* http://www.vibrationdata.com/math/Laplace_Transforms.pdf */
554 /* combined shifting and scaling */
556 (assuming (v - a*b > 0,
557 laplace (exp (a * b * u) * foo (a * u), u, v)));
558 1/a * 'at ('laplace (foo (u), u, v), v = v/a - b);
560 laplace (exp (a * b * u) * cosh (a * u), u, v);
561 ''(ev (1/a * 'at ('laplace (foo (u), u, v), v = v/a - b), 'laplace (foo (u), u, v) = v/(v^2 - 1), nouns, ratsimp));
565 laplace ('diff (delta (u), u), u, v);
573 makefact (laplace (u^(n - 1)/(n - 1)!, u, v));
576 makefact (laplace (u^n, u, v));
579 laplace (u^(a - 1) / gamma (a), u, v);
582 (mycompare (actual, expected) := if ratsimp (actual - expected) = 0 then true else actual, 0);
586 (laplace (u * exp (-a * u), u, v),
591 (makefact (laplace (u^(n - 1)/(n - 1)! * exp (-a * u), u, v)),
596 (laplace (1 - exp (-a * u), u, v),
601 (laplace (1/(b - a) * (exp (-a * u) - exp (-b * u)), u, v),
602 1/((v + a)*(v + b)));
606 (laplace (1/(a*b) + exp (-a * u)/(a * (a - b)) + exp (-b * u)/(b * (b - a)), u, v),
607 1/(v * (v + a) * (v + b)));
611 (laplace (1/(a - b) * (a * exp (-a * u) - b * exp (-b * u)), u, v),
612 v/((v + a) * (v + b)));
616 (laplace ((1 - a*u) * exp (-a * u), u, v),
621 (laplace (1/a^2 * (1 - (1 + a*u) * exp (-a * u)), u, v),
626 (laplace (sin (a * u + b), u, v),
627 (sin (b) * v + cos (b) * a) / (v^2 + a^2));
631 (laplace (u * cos (a*u), u, v),
632 (v^2 - a^2)/(v^2 + a^2)^2);
636 (laplace (1/a^2 * (1 - cos (a*u)), u, v),
637 1/(v * (v^2 + a^2)));
640 /* resume here with eq 2.24 */
642 /* http://www.stanford.edu/~boyd/ee102/laplace-table.pdf */
646 /* bug reported to mailing list 2012-11-25: laplace prevents user-defined simplification */
648 (kill (foo), tellsimpafter ('laplace (foo (w), w, x), FOO (x)), laplace (foo (w), w, x));
651 /* laplace(t^z,t,s) did not properly support complex z expressions
653 * See mailing list thread "Laplace transform for complex exponent polynomial functions."
657 laplace (t^(v*%i), t, s);
658 s^((-%i*v)-1)*gamma(%i*v+1);
660 assuming (u > -1, laplace (t^(u + v*%i), t, s));
661 s^((-%i*v)-u-1)*gamma(%i*v+u+1);
664 assuming(a > 0, laplace(delta(t - a) * sin(b*t), t, s));
665 -(((%i*%e^(2*%i*a*b)-%i)*%e^(-(a*s)-%i*a*b))/2);
668 declare (z, complex), assume (realpart (z) > -1), 0);
672 s^((-z)-1)*gamma(z+1);
674 (remove (z, complex), forget (realpart (z) > -1), 0);
678 (remove (n, integer), killcontext (mycontext));
681 (reset(assume_pos_pred),1);
684 (reset(assume_pos),1);