MaplePrimes Posts

MaplePrimes Posts are for sharing your experiences, techniques and opinions about Maple, MapleSim and related products, as well as general interests in math and computing.

Latest Post
  • Latest Posts Feed
  • I have noticed that there exists a Stack Exchange site for mathematica, and not for maple. My discussions with the part of Stack Exchange that handle the creation of a new Stack Exchange community have said that I must accrue a certain level of interest in the subject in order for it to be approved, and so I thought I would begin here to see if there is suffice level of interest.

    This would not diminish the use of the Maple Primes forum, and an additional proposal, in consideration of the years od dedication that have gone into this domain, be to pool the data between the two, make reputation points the same on both, perhaps even user profiles and questions answered already linkable, and all of the questions already addressed here showing up in the search on both domains.

    I am proposing this simply because I want to encourage the use of maple, and have noted that Stack Exchange is very popular. 

    So I am posting this to get overall feedback from other Maple users, as to what their opinion is regarding this proposal, and advice on whether it should and how it ought to be pursued.

    Integral Transforms (revamped) and PDEs

     

    Integral transforms, implemented in Maple as the inttrans  package, are special integrals that appear frequently in mathematical-physics and that have remarkable properties. One of the main uses of integral transforms is for the computation of exact solutions to ordinary and partial differential equations with initial/boundary conditions. In Maple, that functionality is implemented in dsolve/inttrans  and in pdsolve/boundary conditions .

     

    During the last months, we have been working heavily on several aspects of these integral transform functions and this post is about that. This is work in progress, in collaboration with Katherina von Bulow

     

    The integral transforms are represented by the commands of the inttrans  package:

    with(inttrans)

    [addtable, fourier, fouriercos, fouriersin, hankel, hilbert, invfourier, invhilbert, invlaplace, invmellin, laplace, mellin, savetable, setup]

    (1)

    Three of these commands, addtable, savetable, and setup (this one is new, only present after installing the Physics Updates) are "administrative" commands while the others are computational representations for integrals. For example,

    FunctionAdvisor(integral_form, fourier)

    [fourier(a, b, z) = Int(a/exp(I*b*z), b = -infinity .. infinity), MathematicalFunctions:-`with no restrictions on `(a, b, z)]

    (2)

    FunctionAdvisor(integral_form, mellin)

    [mellin(a, b, z) = Int(a*b^(z-1), b = 0 .. infinity), MathematicalFunctions:-`with no restrictions on `(a, b, z)]

    (3)

    For all the integral transform commands, the first argument is the integrand, the second one is the dummy integration variable of a definite integral and the third one is the evaluation point. (also called transform variable). The integral representation is also visible using the convert network

    laplace(f(t), t, s); % = convert(%, Int)

    laplace(f(t), t, s) = Int(f(t)*exp(-s*t), t = 0 .. infinity)

    (4)

    Having in mind the applications of these integral transforms to compute integrals and exact solutions to PDE with boundary conditions, five different aspects of these transforms received further development:

    • 

    Compute Derivatives: Yes or No

    • 

    Numerical Evaluation

    • 

    Two Hankel Transform Definitions

    • 

    More integral transform results

    • 

    Mellin and Hankel transform solutions for Partial Differential Equations with boundary conditions


    The project includes having all these tranforms available at user level (not ready), say as FourierTransform for inttrans:-fourier, so that we don't need to input with(inttrans) anymore. Related to these changes we also intend to have Heaviside(0) not return undefined anymore, and return itself instead, unevaluated, so that one can set its value according to the problem/preferred convention (typically 0, 1/2 or 1) and have all the Maple library following that choice.

    The material presented in the following sections is reproducible already in Maple 2019 by installing the latest Physics Updates (v.435 or higher),

    Compute derivatives: Yes or No.

     

    For historical reasons, previous implementations of these integral transform commands did not follow a standard paradigm of computer algebra: "Given a function f(x), the input diff(f(x), x) should return the derivative of f(x)". The implementation instead worked in the opposite direction: if you were to input the result of the derivative, you would receive the derivative representation. For example, to the input laplace(-t*f(t), t, s) you would receive d*laplace(f(t), t, s)/ds. This is particularly useful for the purpose of using integral transforms to solve differential equations but it is counter-intuitive and misleading; Maple knows the differentiation rule of these functions, but that rule was not evident anywhere. It was not clear how to compute the derivative (unless you knew the result in advance).

     

    To solve this issue, a new command, setup, has been added to the package, so that you can set "whether or not" to compute derivatives, and the default has been changed to computederivatives = true while the old behavior is obtained only if you input setup(computederivatives = false). For example, after having installed the Physics Updates,

    Physics:-Version()

    `The "Physics Updates" version in the MapleCloud is 435 and is the same as the version installed in this computer, created 2019, October 1, 12:46 hours, found in the directory /Users/ecterrab/maple/toolbox/2019/Physics Updates/lib/`

    (1.1)

    the current settings can be queried via

    setup(computederivatives)

    computederivatives = true

    (1.2)

    and so differentiating returns the derivative computed

    (%diff = diff)(laplace(f(t), t, s), s)

    %diff(laplace(f(t), t, s), s) = -laplace(f(t)*t, t, s)

    (1.3)

    while changing this setting to work as in previous releases you have this computation reversed: you input the output (1.3) and you get the corresponding input

    setup(computederivatives = false)

    computederivatives = false

    (1.4)

    %diff(laplace(f(t), t, s), s) = -laplace(t*f(t), t, s)

    %diff(laplace(f(t), t, s), s) = diff(laplace(f(t), t, s), s)

    (1.5)

    Reset the value of computederivatives

    setup(computederivatives = true)

    computederivatives = true

    (1.6)

    %diff(laplace(f(t), t, s), s) = -laplace(t*f(t), t, s)

    %diff(laplace(f(t), t, s), s) = -laplace(f(t)*t, t, s)

    (1.7)

    In summary: by default, derivatives of integral transforms are now computed; if you need to work with these derivatives as in  previous releases, you can input setup(computederivatives = false). This setting can be changed any time you want within one and the same Maple session, and changing it does not have any impact on the performance of intsolve, dsolve and pdsolve to solve differential equations using integral transforms.

    ``

    Numerical Evaluation

     

    In previous releases, integral transforms had no numerical evaluation implemented. This is in the process of changing. So, for example, to numerically evaluate the inverse laplace transform ( invlaplace  command), three different algorithms have been implemented: Gaver-Stehfest, Talbot and Euler, following the presentation by Abate and Whitt, "Unified Framework for Numerically Inverting Laplace Transforms", INFORMS Journal on Computing 18(4), pp. 408–421, 2006.

     

    For example, consider the exact solution to this partial differential equation subject to initial and boundary conditions

    pde := diff(u(x, t), x) = 4*(diff(u(x, t), t, t))

    iv := u(x, 0) = 0, u(0, t) = 1

     

    Note that these two conditions are not entirely compatible: the solution returned cannot be valid for x = 0 and t = 0 simultaneously. However, a solution discarding that point does exist and is given by

    sol := pdsolve([pde, iv])

    u(x, t) = -invlaplace(exp(-(1/2)*s^(1/2)*t)/s, s, x)+1

    (2.1)

    Verifying the solution, one condition remains to be tested

    pdetest(sol, [pde, iv])

    [0, 0, -invlaplace(exp(-(1/2)*s^(1/2)*t)/s, s, 0)]

    (2.2)

    Since we now have numerical evaluation rules, we can test that what looks different from 0 in the above is actually 0.

    zero := [0, 0, -invlaplace(exp(-(1/2)*s^(1/2)*t)/s, s, 0)][-1]

    -invlaplace(exp(-(1/2)*s^(1/2)*t)/s, s, 0)

    (2.3)

    Add a small number to the initial value of t to skip the point t = 0

    plot(zero, t = 0+10^(-10) .. 1)

     

    The default method used is the method of Euler sums and the numerical evaluation is performed as usual using the evalf command. For example, consider

    F := sin(sqrt(2*t))

     

    The Laplace transform of F is given by

    LT := laplace(F, t, s)

    (1/2)*2^(1/2)*Pi^(1/2)*exp(-(1/2)/s)/s^(3/2)

    (2.4)

    and the inverse Laplace transform of LT in inert form is

    ILT := %invlaplace(LT, s, t)

    %invlaplace((1/2)*2^(1/2)*Pi^(1/2)*exp(-(1/2)/s)/s^(3/2), s, t)

    (2.5)

    At t = 1 we have

    eval(ILT, t = 1)