MaplePrimes - Questions, Posts, help, musings, answers ... all things Maple and math

Ask a Question

Create a Post

Featured Post

ecterrab

7596

Integral Transforms (revamped...

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:

(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,

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

(2)

[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) = 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 anymore. Related to these changes we also intend to have not return 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 , the input should return the derivative of ". 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 you would receive . 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 . For example, after having installed the Physics Updates,

`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

(1.2)

and so differentiating returns the derivative computed

(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

(1.4)

(1.5)

Reset the value of computederivatives

(1.6)

(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 . 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

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

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

[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.

-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

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

The Laplace transform of F is given by

(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

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

(2.5)

At we have

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

(2.6)

(2.7)

This result is consistent with the one we get if we first compute the exact form of the inverse Laplace transform at t = 1:

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

(2.8)

(2.9)

In addition to the standard use of evalf to numerically evaluate inverse Laplace transforms, one can invoke each of the three different methods implemented using the MathematicalFunctions:-Evalf command

(2.10)

(2.11)

(2.12)

(2.13)

Regarding the method we use by default: from a numerical experiment with varied problems we have concluded that our implementation of the Euler (sums) method is faster and more accurate than the other two.

Two Hankel transform definitions

In previous Maple releases, the definition of the Hankel transform was given by

hankel(f(t), t, s, nu) = Int(f(t)*sqrt(s*t)*BesselJ(nu, s*t), t = 0 .. infinity)

where is the function. This definition, sometimes called alternative definition of the Hankel transform, has the inconvenience of the square root in the integrand, complicating the form of the hankel transform for the Laplacian in cylindrical coordinates. On the other hand, the definition more frequently used in the literature is

hankel(f(t), t, s, nu) = Int(f(t)*t*BesselJ(nu, s*t), t = 0 .. infinity)

With it, the Hankel transform of diff(u(r, t), r, r)+(diff(u(r, t), r))/r+diff(u(r, t), t, t) is given by the simple ODE form d^2*`&Hopf;`(k, t)/dt^2-k^2*`&Hopf;`(k, t) . Not just this transform but several other ones acquire a simpler form with the definition that does not have a square root in the integrand.

So the idea is to align Maple with this simpler definition, while keeping the previous definition as an alternative. Hence, by default, when you load the inttrans package, the new definition in use for the Hankel transform is

hankel(f(t), t, s, nu) = Int(f(t)*t*BesselJ(nu, s*t), t = 0 .. infinity)

(3.1)

You can change this default so that Maple works with the alternative definition as in previous releases. For that purpose, use the new inttrans:-setup command (which you can also use to query about the definition in use at any moment):

(3.2)

This change in definition is automatically taken into account by other parts of the Maple library using the Hankel transform. For example, the differentiation rule with the new definition is

(3.3)

This differentiation rule resembles (is connected to) the differentiation rule for BesselJ, and this is another advantage of the new definition.

%diff(BesselJ(nu, z), z) = -BesselJ(nu+1, z)+nu*BesselJ(nu, z)/z

(3.4)

Furthermore, several transforms have acquired a simpler form, as for example:

%hankel(exp(I*a*r)/r, r, k, 0) = 1/(-a^2+k^2)^(1/2)

(3.5)

Let's compare: make the definition be as in previous releases.

(3.6)

hankel(f(t), t, s, nu) = Int(f(t)*(s*t)^(1/2)*BesselJ(nu, s*t), t = 0 .. infinity)

(3.7)

The differentiation rule with the previous (alternative) definition was not as simple:

(3.8)

And the transform (3.5) was also not so simple:

%hankel(exp(I*a*r)/r, r, k, 0) = (I*a*hypergeom([3/4, 3/4], [3/2], a^2/k^2)*GAMMA(3/4)^4+Pi^2*k*hypergeom([1/4, 1/4], [1/2], a^2/k^2))/(k*Pi*GAMMA(3/4)^2)

(3.9)

Reset to the new default value of the definition.

(3.10)

hankel(f(t), t, s, nu) = Int(f(t)*t*BesselJ(nu, s*t), t = 0 .. infinity)

(3.11)

More integral transform results

The revision of the integral transforms includes also filling gaps: previous transforms that were not being computed are now computed. Still with the Hankel transform, consider the operators

`D/t` := proc (u) options operator, arrow; (diff(u, t))/t end proc

Being able to transform these operators into algebraic expressions or differential equations of lower order is key for solving PDE problems with Boundary Conditions.

Consider, for instance, this ODE

(4.1)

((diff(diff(diff(u(t), t), t), t))*t^3-12*(diff(diff(u(t), t), t))*t^2+57*(diff(u(t), t))*t-105*u(t))/t^3

(4.2)

Its Hankel transform is a simple algebraic expression

(4.3)

An example with formula_plus

((diff(diff(diff(diff(u(t), t), t), t), t))*t^4+38*(diff(diff(diff(u(t), t), t), t))*t^3+477*(diff(diff(u(t), t), t))*t^2+2295*(diff(u(t), t))*t+3465*u(t))/t^4

(4.4)

(4.5)

In the case of hankel , not just differential operators but also several new transforms are now computable

piecewise(nu = 0, Dirac(k)/k, nu/k^2)

(4.6)

piecewise(And(nu = 0, m = 0), Dirac(k)/k, 2^(m+1)*k^(-m-2)*GAMMA(1+(1/2)*m+(1/2)*nu)/GAMMA((1/2)*nu-(1/2)*m))

(4.7)

Mellin and Hankel transform solutions for Partial Differential Equations with Boundary Conditions

In previous Maple releases, the Fourier and Laplace transforms were used to compute exact solutions to PDE problems with boundary conditions. Now, Mellin and Hankel transforms are also used for that same purpose.

Example:

pde := x^2*(diff(u(x, y), x, x))+x*(diff(u(x, y), x))+diff(u(x, y), y, y) = 0

iv := u(x, 0) = 0, u(x, 1) = piecewise(0 <= x and x < 1, 1, 1 < x, 0)

(5.1)

As usual, you can let pdsolve choose the solving method, or indicate the method yourself:

pde := diff(u(r, t), r, r)+(diff(u(r, t), r))/r+diff(u(r, t), t, t) = -Q__0*q(r)
iv := u(r, 0) = 0

u(r, t) = Q__0*(-hankel(exp(-s*t)*hankel(q(r), r, s, 0)/s^2, s, r, 0)+hankel(hankel(q(r), r, s, 0)/s^2, s, r, 0))

(5.2)

It is sometimes preferable to see these solutions in terms of more familiar integrals. For that purpose, use

u(r, t) = Q__0*(-(Int(exp(-s*t)*(Int(q(r)*r*BesselJ(0, r*s), r = 0 .. infinity))*BesselJ(0, r*s)/s, s = 0 .. infinity))+Int((Int(q(r)*r*BesselJ(0, r*s), r = 0 .. infinity))*BesselJ(0, r*s)/s, s = 0 .. infinity))

(5.3)

An example where the hankel transform is computable to the end:

pde := c^2*(diff(u(r, t), r, r)+(diff(u(r, t), r))/r) = diff(u(r, t), t, t)
iv := u(r, 0) = A*a/(a^2+r^2)^(1/2), (D[2](u))(r, 0) = 0

u(r, t) = (1/2)*A*a*((-c^2*t^2+(2*I)*a*c*t+a^2+r^2)^(1/2)+(-c^2*t^2-(2*I)*a*c*t+a^2+r^2)^(1/2))/((-c^2*t^2-(2*I)*a*c*t+a^2+r^2)^(1/2)*(-c^2*t^2+(2*I)*a*c*t+a^2+r^2)^(1/2))

(5.4)

Download Integral_Transforms_(revamped).mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

October 01

9 14

Read Full Post

Featured Post

ecterrab

7596

Splitting PDE parameterized...

Splitting PDE parameterized symmetries

and Parameter-continuous symmetry transformations

The determination of symmetries for partial differential equation systems (PDE) is relevant in several contexts, the most obvious of which is of course the determination of the PDE solutions. For instance, generally speaking, the knowledge of a N-dimensional Lie symmetry group can be used to reduce the number of independent variables of PDE by N. So if PDE depends only on N independent variables, that amounts to completely solving it. If only N-1 symmetries are known or can be successfully used then PDE becomes and ODE; etc., all advantageous situations. In Maple, a complete set of symmetry commands, to perform each step of the symmetry approach or several of them in one go, is part of the PDEtools package.

Besides the dependent and independent variables, PDE frequently depends on some constant parameters, and besides the PDE symmetries for arbitrary values of those parameters, for some particular values of them, PDE transforms into a completely different problem, admitting different symmetries. The question then is: how can you determine those particular values of the parameters and the corresponding different symmetries? That was the underlying subject of a recent question in Mapleprimes. The answer to those questions is relatively simple and yet not entirely obvious for most of us, motivating this post, organized briefly around one example.

To reproduce the input/output below you need Maple 2019 and to have installed the Physics Updates v.449 or higher.

Consider the family of Korteweg-de Vries equation for involving three constant parameters . For convenience (simpler input and more readable output) use the diff_table and declare commands

(1)

(2)

This pde admits a 4-dimensional symmetry group, whose infinitesimals - for arbitrary values of the parameters - are given by

(3)

Looking at pde (1) as a nonlinear problem in u, a, b and q, it splits into four cases for some particular values of the parameter:

>	$pde__cases := casesplit(bu(x, t)(diff(u(x, t), x))+a(diff(u(x, t), x))+q(diff(diff(diff(u(x, t), x), x), x))+diff(u(x, t), t) = 0, parameters = {a, b, q}, caseplot)$

`casesplit/ans`([diff(diff(diff(u(x, t), x), x), x) = -(b*u(x, t)*(diff(u(x, t), x))+a*(diff(u(x, t), x))+diff(u(x, t), t))/q], [q <> 0]), `casesplit/ans`([diff(u(x, t), x) = -(diff(u(x, t), t))/(b*u(x, t)+a), q = 0], [b*u(x, t)+a <> 0]), `casesplit/ans`([u(x, t) = -a/b, q = 0], [b <> 0]), `casesplit/ans`([diff(u(x, t), t) = 0, a = 0, b = 0, q = 0], [])

(4)

The legend above indicates the pivots and the tree of cases, depending on whether each pivot is equal or different from 0. At the end there is the algebraic sequence of cases. The first case is the general case, for which the symmetry infinitesimals were computed as above, but clearly the other three cases admit more general symmetries. Consider for instance the second case, pass the to ignore the parameterizing equation , and you get

[_xi[x](x, t, u) = _F3(x, t, u), _xi[t](x, t, u) = Intat(((b*u+a)*(D[1](_F3))(_a, ((b*u+a)*t-x+_a)/(b*u+a), u)-_F1(u, ((b*u+a)*t-x)/(b*u+a))*b+(D[2](_F3))(_a, ((b*u+a)*t-x+_a)/(b*u+a), u))/(b*u+a)^2, _a = x)+_F2(u, ((b*u+a)*t-x)/(b*u+a)), _eta[u](x, t, u) = _F1(u, ((b*u+a)*t-x)/(b*u+a))]

(5)

These infinitesimals are indeed much more general than , in fact so general that (5) is almost unreadable ... Specialize the three arbitrary functions into something "easy" just to be able follow - e.g. take _F1 to be just the + operator, _F2 the * operator and

[_xi[x](x, t, u) = 1, _xi[t](x, t, u) = Intat(-(u+((b*u+a)*t-x)/(b*u+a))*b/(b*u+a)^2, _a = x)+u*((b*u+a)*t-x)/(b*u+a), _eta[u](x, t, u) = u+((b*u+a)*t-x)/(b*u+a)]

(6)

[_xi[x](x, t, u) = 1, _xi[t](x, t, u) = (b^3*t*u^4+((3*a*t-x)*u^3-u^2*x-t*x*u)*b^2+((3*a^2*t-2*a*x)*u^2-a*u*x-a*t*x+x^2)*b+a^2*u*(a*t-x))/(b*u+a)^3, _eta[u](x, t, u) = (b*u^2+(b*t+a)*u+a*t-x)/(b*u+a)]

(7)

This symmetry is of course completely different than computed for the general case.

The symmetry (7) can be verified against or directly against pde after substituting .

(8)

(9)

(10)

Summarizing: "to split PDE symmetries into cases according to the values of the PDE parameters, split the PDE into cases with respect to these parameters (command PDEtools:-casesplit ) then compute the symmetries for each case"

Parameter continuous symmetry transformations

A different, however closely related question, is whether pde admits "symmetries with respect to the parameters a, b and q", so whether exists continuous transformations of the parameters a, b and q that leave pde invariant in form.

Beforehand, note that since the parameters are constants with regards to the dependent and independent variables (here ), such continuous symmetry transformations cannot be used directly to compute a solution for pde. They can, however, be used to reduce the number of parameters. And in some contexts, that is exactly what we need, for example to entirely remove the splitting into cases due to their presence, or to proceed applying a solving method that is valid only when there are no parameters (frequently the case when computing exact solutions to "PDE & Boundary Conditions").

To compute such "continuous symmetry transformations of the parameters" that leave pde invariant one can always think of these parameters as "additional independent variables of pde". In terms of formulation, that amounts to replacing the dependency in the dependent variable, i.e. replace by

(11)

Compute now the infinitesimals: note there are now three additional ones, related to continuous transformations of and - for readability, avoid displaying the redundant functionality on the left-hand sides of these infinitesimals

[_xi[x] = (1/3)*(_F4(a, b, q)*q+_F3(a, b, q))*x/q+_F6(a, b, q)*t+_F7(a, b, q), _xi[t] = _F4(a, b, q)*t+_F5(a, b, q), _xi[a] = _F1(a, b, q), _xi[b] = _F2(a, b, q), _xi[q] = _F3(a, b, q), _eta[u] = (1/3)*((b*u+a)*_F3(a, b, q)-2*((b*u+a)*_F4(a, b, q)+(3/2)*u*_F2(a, b, q)+(3/2)*_F1(a, b, q)-(3/2)*_F6(a, b, q))*q)/(b*q)]

(12)

This result is more general than what is convenient for algebraic manipulations, so specialize the seven arbitrary functions of and keep only the first symmetry that result from this specialization: that suffices to illustrate the removal of any of the three parameters a, b, or q

(13)

To remove the parameters, as it is standard in the symmetry approach, compute a transformation to canonical coordinates, with respect to the parameter a. That means a transformation that changes the list of infinitesimals, or likewise its infinitesimal generator representation,

proc (f) options operator, arrow; diff(f, a)-(diff(f, u))/b end proc

(14)

into or its equivalent generator representation

That same transformation, when applied to , entirely removes the parameter a.

The transformation is computed using CanonicalCoordinates and the last argument indicates the "independent variable" (in our case a parameter) that the transformation should remove. We choose to remove a

(15)

(16)

Invert this transformation in order to apply it

(17)

The next step is not necessary, but just to understand how all this works, verify its action over the infinitesimal generator proc (f) options operator, arrow; diff(f, a)-(diff(f, u))/b end proc

ChangeSymmetry({a = alpha, b = beta, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*beta-alpha)/beta}, proc (f) options operator, arrow; diff(f, a)-(diff(f, u))/b end proc, [upsilon(xi, tau, alpha, beta, chi), xi, tau, alpha, beta, chi])

(18)

Now that we see the transformation (17) is the one we want, just use it to change variables in

>	$PDEtools:-dchange({a = alpha, b = beta, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*beta-alpha)/beta}, pde__xtabq, [upsilon(xi, tau, alpha, beta, chi), xi, tau, alpha, beta, chi], simplify)$

(19)

As expected, this result depends only on two parameters, , and , and the one equivalent to a (that is , see the transformation used (17)), is not present anymore.

To remove b or q we use the same steps, (15), (17) and (19), just changing the parameter to be removed, indicated as the last argument in the call to CanonicalCoordinates . For example, to eliminate b (represented in the new variables by ), input

(20)

(21)

>	$PDEtools:-dchange({a = beta, b = alpha, q = chi, t = tau, x = xi, u(x, t, a, b, q) = (upsilon(xi, tau, alpha, beta, chi)*alpha-beta)/alpha}, pde__xtabq, [upsilon(xi, tau, alpha, beta, chi), xi, tau, alpha, beta, chi], simplify)$

(22)

and as expected this result does not contain To remove a second parameter, the whole cycle is repeated starting with computing infinitesimals, for instance for (22). Finally, the case of function parameters is treated analogously, by considering the function parameters as additional dependent variables instead of independent ones.

Download How_to_split_symmetries_into_cases_(II).mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

October 18

5 0

Read Full Post

Maple May benefit from starting a Stack Exchange...

by: Adam Ledger 280 Product: Maple

October 12 2019

4 7

Pole simulation with tensions in MapleSim

by: Lenin Araujo Castillo 1480 Products: Maple 2018 MapleSim 2018

October 14 2019

2 0

MaplePrimes Questions

Recent

Unanswered

Maple

MapleSim

Share via:

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

Featured Post

Integral Transforms (revamped...

Featured Post

Splitting PDE parameterized...

Maple May benefit from starting a Stack Exchange...

Pole simulation with tensions in MapleSim

MaplePrimes Questions

Recent

Unanswered

Maple

MapleSim

Please need help with my proc

parameters assigned in the equation

i can't find our where my error is...

How to improve physics Simplify for an Angular...

How to draw line following the gradient curve

Unexpected difference between products and sums...

Need help with rtable_eval

idea for trigonometric animation in Maple ?

Maple evaluating function when it is impossible

How can I extract DE coefficients to form a system...

display not displaying all plots on one plot

how do i make ( isentropic-process) ?

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again

Ask a Question

Create a Post

Featured Post

Featured Post

MaplePrimes Questions Recent Unanswered Maple MapleSim

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

MaplePrimes Questions

Recent

Unanswered

Maple

MapleSim