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

Ask a Question

Create a Post

MaplePrimes Announcement

Update to Maple 2025 and Announcing MapleSim 2025

by: TechnicalSupport 760 Products: Maple 2025 , MapleSim 2025

2 hours ago

1 0

Maple 2025.1

We have just released an update to Maple. Maple 2025.1 includes several enhancements to the new interface, as well as various small corrections throughout the product. As always, we recommend that all Maple 2025 users install this update.

In particular, please note that this update includes a fix to the problem where new documents were opening in a new window instead of a new tab. Thanks for helping us, and other users, by letting us know!

This update is available through Tools>Check for Updates in Maple, and is also available from the Maple 2025.1 download page on web site, where you can also find more details.

MapleSim 2025

We are happy to announce that we just released MapleSim 2025. This release includes a new component library to support the modeling of motor drives and updates to several in-product apps that make it even easier to perform optimization and analysis.

See What’s New in MapleSim for details.

Read Full Post

Featured Post

dharr

8070

Laplace Adomian Decomposition method for solving PDEs

Thanks to @salim-barzani, I became interested in the Laplace Adomian Decompositon method to solve partial differential equations. It produces a sum of analytical components that converge to the exact solution. Probably it is not that efficient for numerical solutions, but the possibility of finding an exact solution to nonlinear PDEs by finding a formula for the components and therefore the infinite sum is interesting.

The version here covers many of the common forms of PDEs, but is still a work in progress. The method for finding the Adomian polynomials could be improved by using one of the reccurence formulas in the literature. Extension to more PDE forms, ODEs etc are possible and I will attempt some of these before uploading an improved version to the application centre. In the meantime, feedback about bugs, usability etc. are appreciated.

Laplace Adomian Decomposition method to solve partial differential equations.
D.A. Harrington, June 2025, v 1.16.
Procedure LAD is in the startup code region.

Arguments:

Partial differential equation in one "time" variable and several other variables (called spatial variables here). Specified as an equation, or an expression implicitly equal to zero.
Comprising: (1) one time derivative term of order (denoted ), (2) linear terms that are derivatives in the spatial variables (), (3) nonlinear terms with derivatives in the spatial variables (). Together, (1)-(3) is a multivariate polynomial in the dependent variable, its derivatives, and any parameters. (4) inhomogenous terms in the time and spatial variables that do not contain the dependent variable or its derivatives (). The general form is .
Any symbolic parameters are assumed to be real. Any numeric coefficients or parameters should be of type algebraic. Floating point values will be converted to rationals with convert(value, rational).

2.	Function to solve for, e.g., .

3.	initial condition(s) (at time zero). For , these must be in a list, and are the values of the function and its successive derivatives with respect to time, evaluated at time zero.

4.	name of time variable.

5.	number of iterations.

(optional) order for a fractional derivative, specified as $fracorder = name$ or $fracorder = numeric*value$ (floating point values are converted to type rational). The permissible values are related to the value of : , i.e., an th order time derivative is specified in the pde and initial conditions are provided, as well as the value of in the correct range. A symbolic is assumed to be in this range.

infolevel[LAD] may be set to different values to print out additional information as follows. Greater numbers include the information from smaller values.

1.	The nonlinear expansion variables (those in the Adomian polynomials).

2.	The different parts of the the pde () in jet notation.

3.	Progress is indicated, as the time at each iteration.

4.	Values of the components of the solution as they are produced (one per iteration).

It may be useful to load the Physics package if derivatives contain abs. LAD converts these to a form without abs but with the conjugate of the dependent variable. The Physics package affects the internal simplifications and may affect the form of the output or the running time.

Example from A-M. Wazwaz, Applied Mathematics and Computation 111 (2000) 53. Sec. 4.

LAD: L = U[t]; R = 0; N = -U^2*U[x]; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: [U, U[x]]

By extrapolating to an infinite number of terms, we can find an exact solution. (Wazwaz calculates the coefficient of incorrectly and deduces an incorrect exact solution.) Multiplying by gives a series in for which guessgf can use the coefficients to find the exact solution.

algsubs(x*t = y, expand(approx*t)); [seq(coeff(%, y, i), i = 0 .. degree(%, y))]; gfun:-guessgf(%, y)[1]; exact := (eval(%, y = x*t))/t

Check

An example with inhomogeneous terms (but no nonlinear part). Ex. 3 from R. Shah et al, Entropy 21 (2019) 335.

pde := diff(U(x, t), t)+diff(U(x, t), `$`(x, 3)) = -sin(Pi*x)*sin(t)-Pi^3*cos(Pi*x)*cos(t); inx := sin(Pi*x); exact := sin(Pi*x)*cos(t); pdetest(U(x, t) = exact, [pde, U(x, 0) = inx])

LAD: L = U[t]; R = -U[x,x,x]; N = 0; G = -sin(Pi*x)*sin(t)-Pi^3*cos(Pi*x)*cos(t).

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: []

Expand the exact solution as a series in to the same order and verify that it is the same as above.

An example with a complex solution

pde := I*(diff(U(x, t), t))+diff(U(x, t), `$`(x, 2))+2*abs(U(x, t))^2*U(x, t); inx := sech(x); exact := sech(x)*exp(I*t); `assuming`([simplify(pdetest(U(x, t) = exact, [pde, U(x, 0) = inx]))], [real])

LAD: L = I*U[t]; R = -U[x,x]; N = -2*U^2*C; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: [C, U]

sech(x)+I*sech(x)*t-(1/2)*sech(x)*t^2-((1/6)*I)*sech(x)*t^3+(1/24)*sech(x)*t^4+((1/120)*I)*sech(x)*t^5-(1/720)*sech(x)*t^6-((1/5040)*I)*sech(x)*t^7+(1/40320)*sech(x)*t^8+((1/362880)*I)*sech(x)*t^9-(1/3628800)*sech(x)*t^10

The series for is apparent.

An example with fractional order differentiation wrt of order . For , the derivative must be entered with order and there must be a list of initial conditions:

Ex. 1 from R. Shah et al, Entropy 21 (2019) 335

Symbolic order is accepted.

>	$approx := LAD(pde, U(x, t), inx, t, 6, fracorder = alpha)$

LAD: L = U[t]; R = -2*U[x]-U[x,x,x]; N = 0; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: []

sin(x)-cos(x)*t^alpha/GAMMA(1+alpha)-sin(x)*t^(2*alpha)/GAMMA(1+2*alpha)+cos(x)*t^(3*alpha)/GAMMA(1+3*alpha)+sin(x)*t^(4*alpha)/GAMMA(1+4*alpha)-cos(x)*t^(5*alpha)/GAMMA(1+5*alpha)-sin(x)*t^(6*alpha)/GAMMA(1+6*alpha)

Giving a numerical value for the order will generally be more efficient. For , we can find an exact solution.

>	$approx := collect(LAD(pde, U(x, t), inx, t, 20, fracorder = 1/2), [sin, cos])$

LAD: L = U[t]; R = -2*U[x]-U[x,x,x]; N = 0; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: []

The first series is , the second series is not so obvious, After some scaling, the second series may be passed to guessgf

ser2 := `assuming`([simplify(expand(sqrt(Pi)*op(2, approx)/(sqrt(t)*cos(x))))], [t > 0]); [seq(coeff(ser2, t, i), i = 0 .. 9)]; ser2ex := gfun:-guessgf(%, t)[1]

-2+(4/3)*t-(8/15)*t^2+(16/105)*t^3-(32/945)*t^4+(64/10395)*t^5-(128/135135)*t^6+(256/2027025)*t^7-(512/34459425)*t^8+(1024/654729075)*t^9

-exp(-t)*Pi^(1/2)*erfi(t^(1/2))/t^(1/2)

So the full solution is

Check. fracdiff's default method returns an integral, method = laplace fails, but method = series works with cancellation of all terms but the "left over" last half-integral-order one.

>	$Ord := 20; fracdiff(exact, t, 1/2, method = series, methodoptions = [order = Ord])-series(eval(rhs(pde), U(x, t) = exact), t, Ord)$

-(1048576/319830986772877770815625)*sin(x)*t^(39/2)/Pi^(1/2)+O(t^(39/2))-O(t^20)

A second-order example (linear)

pde := diff(U(x, t), t, t) = -c^2*(diff(U(x, t), `$`(x, 2))); inx := [exp(I*x/c), exp(I*x/c)]; exact := exp(t+I*x/c)

LAD: L = U[t,t]; R = -c^2*U[x,x]; N = 0; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: []

exp(I*x/c)*(1+t)+(1/6)*exp(I*x/c)*t^2*(3+t)+(1/120)*exp(I*x/c)*t^4*(t+5)+(1/5040)*exp(I*x/c)*t^6*(7+t)+(1/362880)*exp(I*x/c)*t^8*(9+t)+(1/39916800)*exp(I*x/c)*t^10*(t+11)+(1/6227020800)*exp(I*x/c)*t^12*(13+t)

Fractional order example

>	$approx := collect(LAD(pde, U(x, t), inx, t, 10, fracorder = 3/2), [exp, t])$

LAD: L = U[t,t]; R = -c^2*U[x,x]; N = 0; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: []

There appear to be two separate series in , with integer and half-integer powers.

t_series := approx/exp(I*x/c); t1, t2 := selectremove(proc (w) options operator, arrow; (degree(w, t))::integer end proc, t_series)

The integer-powers series is the series for with every third term missing: missing powers 2, 5, 8, 11, ...

t1sum := simplify(exp(t)-(sum(t^(3*i+2)/factorial(3*i+2), i = 0 .. infinity))); series(%, t, 17)

(2/3)*exp(t)+(1/3)*(cos((1/2)*3^(1/2)*t)+3^(1/2)*sin((1/2)*3^(1/2)*t))*exp(-(1/2)*t)

series(1+t+(1/6)*t^3+(1/24)*t^4+(1/720)*t^6+(1/5040)*t^7+(1/362880)*t^9+(1/3628800)*t^10+(1/479001600)*t^12+(1/6227020800)*t^13+(1/1307674368000)*t^15+(1/20922789888000)*t^16+O(t^17),t,17)

And for the other one, the following scaled series shows coefficients with powers of two divided by double factorials, every third term missing: missing powers 2,5,8,11,...

expand(sqrt(Pi)*t2/(4*t^(3/2))); all := add(2^i*t^i/doublefactorial(2*(i+1)+1), i = 0 .. 9); missing := add(eval(2^i*t^i/doublefactorial(2*(i+1)+1), i = 3*j+2), j = 0 .. 3)

(8192/6190283353629375)*t^13+(4096/213458046676875)*t^12+(1024/316234143225)*t^10+(512/13749310575)*t^9+(128/34459425)*t^7+(64/2027025)*t^6+(16/10395)*t^4+(8/945)*t^3+(2/15)*t+1/3

1/3+(2/15)*t+(4/105)*t^2+(8/945)*t^3+(16/10395)*t^4+(32/135135)*t^5+(64/2027025)*t^6+(128/34459425)*t^7+(256/654729075)*t^8+(512/13749310575)*t^9

t2sum := 4*t^(3/2)*(sum(2^i*t^i/doublefactorial(2*(i+1)+1), i = 0 .. infinity)-(sum(eval(2^i*t^i/doublefactorial(2*(i+1)+1), i = 3*j+2), j = 0 .. infinity)))/sqrt(Pi)

4*t^(3/2)*((1/4)*(Pi^(1/2)*exp(t)-2*t^(1/2)-Pi^(1/2)*erfc(t^(1/2))*exp(t))/t^(3/2)-(4/105)*t^2*hypergeom([1], [3/2, 11/6, 13/6], (1/27)*t^3))/Pi^(1/2)

Putting it together suggests

(1/3)*exp(I*x/c)*(Pi^(1/2)*exp(-(1/2)*t)*sin((1/2)*3^(1/2)*t)*3^(1/2)+Pi^(1/2)*exp(-(1/2)*t)*cos((1/2)*3^(1/2)*t)+3*Pi^(1/2)*exp(t)*erf(t^(1/2))-(16/35)*t^(7/2)*hypergeom([1], [3/2, 11/6, 13/6], (1/27)*t^3)+2*Pi^(1/2)*exp(t)-6*t^(1/2))/Pi^(1/2)

Check.

>	$Ord := 20; fracdiff(exact, t, 3/2, method = series, methodoptions = [order = Ord+1])-series(eval(rhs(pde), U(x, t) = exact), t, Ord)$

O(t^(39/2))-(1048576/319830986772877770815625)*exp(I*x/c)*t^(39/2)/Pi^(1/2)-O(t^(41/2))

Numerical calculation of a soliton

pde := I*(diff(U(x, t), t))+diff(U(x, t), `$`(x, 2))+U(x, t)^2*conjugate(U(x, t)); exact := sqrt(2*B)*sech(sqrt(B)*(2*k*t-x-x__0))*exp(I*(k*x+(-k^2+B)*t)); `assuming`([simplify(pdetest(U(x, t) = exact, pde))], [real, B > 0]); inx := eval(exact, t = 0)

With numerical parameters it runs faster. Float parameters will be converted to rationals.

>	$params := {B = 2, k = -1/2, x__0 = -2}; inxnum := eval(inx, params); exactnum := eval(exact, params)$

LAD: L = I*U[t]; R = -U[x,x]; N = -U^2*C; G = 0.

LAD: nonlinear expansion variables (conjugated functions denoted by C) are: [C, U]

LAD: calculating component 1 at time 0.

LAD: calculating component 2 at time .221

[Edited for brevity]

LAD: calculating component 27 at time 90.414

LAD: calculating component 28 at time 103.870

LAD: completed at time 105.309

For larger and , the approximation diverges at a point that depends on the number of iterations.

Download Laplace_Adomian_Decomposition_procedure8b.mw

June 17

7 7

Read Full Post

Featured Post

Rouben Rostamian

8456

Optical illusion

In Optical-Illusion--Impossible-Prisms, mcarvalho provides a Maple Learn worksheet to illustrate an intriguing optical illusion. Here I do that illustration in Maple, and in 3D!

The graphics below shows 4 bars. Or is it 3 bars?

Download the worksheet, stare at the graphics for short while to grasp the nature of the optical illusion, and then rotate the 3D graph with the mouse and see the illusion fall apart.

optical-illusion.mw

May 27

6 0

Read Full Post