Pdsolve find exact solutions for systems of partial differential equations (pdes)

pdsolve - find exact solutions for systems of partial differential equations (PDEs)
Calling Sequence
pdsolve(PDE_sys, optional_1, optional_2,...)
Parameters
PDE_sys - system of PDEs; it can contain inequations
optional_i - (optional) arguments that can be given in any order and are described as follows
funcs - set or list with indeterminate functions or names - can also be a ranking
HINT = `+` - indicate that the PDE system is to be solved by separating the variables by sum
HINT = `*` - indicate that the PDE system is to be solved by separating the variables by product
HINT = TWS - indicate that the PDE system is to be solved by computing a traveling wave solution as a power series in tanh(x), where xi represents a linear combination of the independent variables
HINT = TWS(foo) - indicate that the PDE system is to be solved by computing a traveling wave solution as a power series in foo, where foo is the name of a mathematical function
singsol=false - avoid the computation of the singular solutions when the system is nonlinear
mindim=N - avoid the computation of solutions if the dimension of the solution space is less than N
parameters=P - P is a list or set of names or functions which are solving variables with priority lower than any other variable
diffalg - request that the differential elimination step for all PDE systems be performed by using the diffalg package
Description
Given a PDE system, possibly including ODEs, algebraic constraints, and inequations, pdsolve returns an exact solution if: 1) the uncoupling of the system does not require more resources than those available in the computer and 2) the pdsolve routines for solving a PDE subsystem in a single unknown as well as the dsolve routines for solving related ODE systems succeed in solving the PDE and ODE subsystems arising in the uncoupling process. Depending on the complexity of the system, step 2) may not be fully successful. However, step 1) is systematic in that, generally speaking, whenever the system is consistent, the uncoupling is feasible.
The solving process
The solving process consists of two steps.
1. The system is first uncoupled by using differential algebra techniques for polynomial systems. When the PDE system, rational in the unknowns and their derivatives, contains non-polynomial coefficients, the uncoupling is performed by rewriting the system in polynomial form by using a differential extension approach (see PDEtools[dpolyform]).
Due to the intrinsic nature of the differential algebra elimination process, this first step always produces (possibly many) PDE subsystems, such that:
- The equations inside each subsystem satisfy all the integrability conditions;
- The union of the non-singular solutions of each subsystem is equal to the general and singular solutions of the original system;
- Provided that the original system is not subdetermined, one of the PDE subsystems depends on a single unknown.
2. The PDE subsystem for a single unknown is solved, and the unknown removed from all other subsystems. Due to nature of the differential elimination process, after removing this single unknown, there exists a PDE subsystem that depends on a single unknown.
Step 2 is repeated until all the PDE subsystems obtained in step 1 are solved, thus arriving at the solution to the original input PDE system.
When solving each of these PDE subsystems involving a single unknown, two situations may arise: the subsystem consists of only one PDE, or it consists of many PDEs. In the former case, pdsolve solves the problem by using older subroutines. In the latter case, pdsolve uses a set of routines that, at each step, change variables introducing differential invariants as new variables. This approach permits solving the problem when the PDE subsystem is essentially nonlinear.
The output
When successful, pdsolve returns a sequence of solution sets, including the singular solutions of the given PDE system. These singular solutions exist only when the system is nonlinear and can be identified by observing that the dimension of the solution space (basically the number of integration constants appearing) is less than the dimension of the general solution.
When pdsolve is partially successful, the output is a sequence of sets and lists. The sets contain solutions while the lists also contain unsolved differential equations. It is valuable to have these partial solutions because they may be easier to solve. For example, the unsolved equations may be ODEs instead of PDEs, less in number, or of lower differential order.
pdsolve can fail to solve any of the PDE subsystems obtained in step 1 of "The solving process" described above. In such a case, the output consists of a sequences of lists, where each list contains a PDE subsystem. The value of this output is that these PDE subsystems represent the uncoupling (also known as triangularization) of the original system. Depending on the PDE system, this triangularization may represent a relevant step towards solving the problem.
The unknowns of the problem and rankings
By default, when only the PDE system is given, pdsolve considers all the unknown functions present in the system as the unknowns of the problem - herein also called the solving variables. This default can be changed by calling pdsolve with an extra argument, specifying the solving variables either as a set or as a list containing the names of these functions or the functions themselves.
When the solving variables are indicated as a set, the ordering chosen by pdsolve for uncoupling the system is chosen by the PDEtools[casesplit] command, which performs the triangularization . When the solving variables are indicated as a list, for instance, [f, g, h], the uncoupling ordering attempted is: to obtain a PDE subsystem depending on only h, another PDE subsystem depending on only {g, h}, and another one depending on {f, g, h}. If such a triangularization is successful, pdsolve first solves the subsystem for h, then the one for {g,h} by substituting the solution found for h and solving the resulting system for g, and finally the one for {f,g,h} by substituting the solutions found for {h,g} and solving the resulting system for f.
The concept of "solving ordering" for the unknowns of the problem is not different from the concept of "ranking" in differential algebra (for a useful discussion on rankings, see rifsimp[ranking]). In fact, the "solving ordering" accepted by pdsolve can be more elaborate. Consider the ones used to indicate a "ranking" when using the diffalg or DEtools[Rif] packages. For example, by indicating the unknowns of the problem as in [{f, g}, h, {j, k}], one is indicating that pdsolve must decide whether f comes before g (or the opposite), h is ranked third, and then pdsolve must decide whether j comes before k. This flexibility is extremely useful when solving nonlinear PDE systems, where slight changes in the solving ordering can lead to solutions in quite different formats or make an otherwise unsolvable (by pdsolve) problem solvable.
In addition to indicating the solving variables or the solving ordering as explained above, one can use the optional argument parameters=P, where P is a set or list of solving variables with less priority, meaning that they always appear at the end of the list of solving variables described in the previous paragraph. When P is given as a list, the rules for choosing a particular ordering inside it are the ones described in the previous paragraph for the main solving variables.
Other optional arguments
When the PDE system is nonlinear, by default, pdsolve's output contains, in separate solution sets, the singular solutions of the system. This involves both computing the singular cases (when uncoupling the system inside PDEtools[casesplit]) and solving them. In many applications, however, the computation of the singular cases is not relevant. Also their computation is sometimes more time consuming than the computation of the general solution. For these reasons, an optional argument is provided to override this behavior, so that by giving singsol=false these singular solutions are not computed.
When a PDE system is solved, the dimension of the solution space (typically indicated by the number of integration constants present) depends on various factors, not only the "apparent differential order" of the system. In some applications, however, it is the case that the solutions of interest are only those that have at least a minimum number of dimensions. In such a case, to discard solutions with dimensions less than N (some positive integer), use the optional mindim=N argument. For nonlinear PDE systems, this option may also speed up the computational process since there is no time consumed in computing and solving the lower-dimension cases. For linear PDE systems, where there is only one solution, pdsolve returns NULL if N is greater than the dimension of the solution.
It is possible to ask pdsolve to solve a PDE system by trying to separate the variables by sum or product. For this purpose, use the optional argument HINT = `+` or HINT = `*`, respectively.
It is possible to ask pdsolve to solve a PDE system by computing a Traveling Wave Solution, as a power series expansion in tanh (default) or one of the following functions: exp, ln, sin, cos, tan, their multiplicative inverses csc, sec and cot, the corresponding six hyperbolic functions, arcsinh, the twelve elliptic Jacobi functions, the WeierstrassP function and also the identity function x -> x, useful to compute purely polynomial solutions when they exist.
By default, pdsolve performs all differential elimination processes for PDE systems by using the DEtools[Rifsimp] package. This default can be changed by using the optional argument diffalg, anywhere in the calling sequence, which causes the computation to use the diffalg package.
For further mathematical conventions related to the output by pdsolve, see the conventions used in solve, dsolve, and int.
Examples
The computations below were performed on a Pentium II-400 laptop PC running Windows. Maple solved all five examples in 36 seconds.
Example 1.
Compute the point symmetries of an ordinary differential equation (ODE), that is, solve the determining PDE system for the infinitesimals of the symmetry generator. Consider example 11 from Kamke's book:
> withApplyFunction(PDEtools,casesplit,declare)
[casesplit, declare]
> withApplyFunction(DEtools,gensys)
[gensys]
> declareApplyFunction(yApplyFunction(x),prime=x)
y(x) will now be displayed as y
derivatives with respect to x of functions of one variable will now be

displayed with '

> 2
> DifferentialD r n
> ode[11]:=--------------- yApplyFunction(x)+a x yApplyFunction(x) =0
> 2
> DifferentialDx
r n
ode[11] := y'' + a x y = 0
The PDE system satisfied by the symmetries, that is, infinitesimals [xi, eta] of the symmetry generator, of the ODE above is given by
> declareApplyFunction((?,?)ApplyFunction(x,y))
xi(x, y) will now be displayed as xi
eta(x, y) will now be displayed as eta
> sys:=[gensysApplyFunction(ode[11],[?,?]ApplyFunction(x,y))];

for _eq in sys do _eq=0 end do

xi[y, y] = 0
eta[y, y] - 2 xi[x, y] = 0
r n
3 a x y xi[y] + 2 eta[x, y] - xi[x, x] = 0
r n r n
eta a x y n r n r n xi a x r y
------------- + 2 a x y xi[x] - a x y eta[y] + ------------ + eta[x, x] = 0
y x
This is a second order linear PDE system, with two unknowns {xi(x,y), eta(x,y)} and four equations. Its general solution is given by
> sol:=pdsolveApplyFunction(sys)
/ _C1 y (2 + r) \
sol := { eta = - -------------, xi = _C1 x }
\ -1 + n /
Since symmetries are defined up to a multiplicative constant, we can drop _C1. Hence, for arbitrary {a, n, r}, the ODE has only one point symmetry.
Solutions to PDE systems can be tested by using pdetest (the command tests whether the solution satisfies each PDE in the system).
> pdetestApplyFunction(sol,sys)
[0, 0, 0, 0]
Now, an interesting question is: Are there other solutions related to particular values of the parameters n and r?
To answer that question one must solve a more difficult, now nonlinear, PDE system by indicating to pdsolve that n and r are parameters.
Also, in this particular problem, from the form of the ODE y''+a*x^r*y^n = 0, the case n = 1 is of no interest since the ODE would become linear, and so, computing its symmetries is the same as computing its solution. We are interested in a solution for different values of n and r but also for n <> 1. The first step is then to add this inequation to the PDE system.
> sys1:=[opApplyFunction(sys),nne1]
Next we indicate to pdsolve that n and r are parameters of the problem, thus arriving at the desired solutions.
> sol1:=pdsolveApplyFunction(sys1,parameters={n,r})
/
|
sol1 := {eta = y (3 _C1 + _C2 x), n = 2, r = -5, xi = x (_C1 + _C2 x)}, <
|
\

/ (8/7) 2 \

2 \-98 y x a _C1 - 147 y a _C2 x - 6 x _C1/ -20
eta = - -------------------------------------------------, n = 2, r = ---,
343 x a 7

\ /
(8/7) | |

xi = _C1 x + _C2 x >, <
| |
/ \

(6/7)
-147 y x a _C1 - 49 y a _C2 x + 12 _C1 x -15

eta = - ---------------------------------------------, n = 2, r = ---,
343 x a 7

\
(6/7) |

xi = _C1 x + _C2 x >, {eta = -_C1 y (2 + r), n = 2, r = r, xi = _C1 x},
|
/

/ y ((_C1 + _C2 x) r + 2 _C1 + 4 _C2 x)

{ eta = -------------------------------------, n = -r - 3, r = r,
\ 4 + r

\ / _C1 y (2 + r) \

xi = x (_C1 + _C2 x) }, { eta = - -------------, n = n, r = r, xi = _C1 x }
/ \ -1 + n /
> mapApplyFunction(pdetest,[sol1],sys1)
[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0],

[0, 0, 0, 0]]

So there exist particular values of n and r for which the system has additional solutions. The solution set with n and r integers and with xi linear in x is in fact a particular case of the general solution computed previously, but the other solution sets are not.
Example 2.
Compute the general solution of the following (linear) overdetermined system involving two PDEs, three unknown functions, one of which depends on 2 variables and the other two depend on only 1 variable.
> declareApplyFunction(FApplyFunction(r,s),HApplyFunction(r),GApplyFunction(s))
F(r, s) will now be displayed as F
H(r) will now be displayed as H
G(s) will now be displayed as G
> [ / 2 \ 2
> [ |PartialD | PartialD
> sys2:=[uminus0|---------- FApplyFunction(r,s)|+---------- FApplyFunction(r,s)+
> [ | 2 | 2
> [ \PartialDr / PartialDs
>
> DifferentialD DifferentialD
> -------------- HApplyFunction(r)+-------------- GApplyFunction(s)+s=0,
> DifferentialDr DifferentialDs
>
> 2 / 2 \
> PartialD | PartialD |
> ---------- FApplyFunction(r,s)+2 |------------------- FApplyFunction(r,s)|+
> 2 \PartialDs PartialDr /
> PartialDr
>
> 2
> PartialD /DifferentialD \
> ---------- FApplyFunction(r,s) - |-------------- HApplyFunction(r)|+
> 2 \DifferentialDr /
> PartialDs
>
> ]
> DifferentialD ]
> -------------- GApplyFunction(s) - r=0]
> DifferentialDs ]
> ]
sys2 := [-F[r, r] + F[s, s] + H[r] + G[s] + s = 0,

F[r, r] + 2 F[r, s] + F[s, s] - H[r] + G[s] - r = 0]

The solution for the unknowns G, H, is given by
> sol:=pdsolveApplyFunction(sys2)
/ 1 2
sol := { F = _F1(s) + _F2(r) + _F3(s - r) - -- r (r - 3 s),
\ 12

1 2 1 2 \

G = - - s - _F1[s] + _C2, H = _F2[r] - - r + _C1 }
4 4 /
and depends on 3 arbitrary functions, _F1(s), _F2(r), _F3(s-r).
> pdetestApplyFunction(sol,sys2)
[0, 0]
Example 3.
Compute the solution of the following overdetermined first order linear system of 3 PDEs in a single unknown f(x,y,z,t).
> [ /PartialD \ 2 /PartialD
> sys3:=[uminus0y |--------- fApplyFunction(x,y,z,t)|+z |--------- f
> [ \PartialDx / \PartialDz
>
> \ /PartialD \ 2
> ApplyFunction(x,y,z,t)|+3 t z |--------- fApplyFunction(x,y,z,t)| - 3 t - 4
> / \PartialDt /
>
> /PartialD \
> fApplyFunction(x,y,z,t) z=0,uminus0y |--------- fApplyFunction(x,y,z,t)| -
> \PartialDy /
>
> /PartialD \ /PartialD
> z |--------- fApplyFunction(x,y,z,t)| - t |--------- fApplyFunction(x,y,z,t)
> \PartialDz / \PartialDt
>
> \ /PartialD \
> |+fApplyFunction(x,y,z,t)=0,uminus0x |--------- fApplyFunction(x,y,z,t)| -
> / \PartialDy /
>
> /PartialD \ ]
> |--------- fApplyFunction(x,y,z,t)|=0];
> \PartialDz / ]

for _eq in sys3 do _eq end do

2 2
-y f[x] + z f[z] + 3 t z f[t] - 3 t - 4 f(x, y, z, t) z = 0
-y f[y] - z f[z] - t f[t] + f(x, y, z, t) = 0
-x f[y] - f[z] = 0
The solving technique makes use of differential invariants and the solution is given by
> sol:=pdsolveApplyFunction(sys3)
/ 2 (3/2) \
| 3 t x t _C1 |
sol := < f(x, y, z, t) = ------- + -------------- >
| x z - y (1/2)|
\ (x z - y) /
> pdetestApplyFunction(sol,sys3)
[0, 0, 0]
Example 4.
Compute the solution of the following nonlinear system, consisting of Burger's equation and a possible potential.
> declareApplyFunction(uApplyFunction(x,t),vApplyFunction(x,t))
u(x, t) will now be displayed as u
v(x, t) will now be displayed as v
> [
> [PartialD /PartialD
> sys4:=[--------- uApplyFunction(x,t)+2 uApplyFunction(x,t) |--------- u
> [PartialDt \PartialDx
> [
>
> / 2 \
> \ |PartialD | PartialD
> ApplyFunction(x,t)| - |---------- uApplyFunction(x,t)|=0,--------- v
> / | 2 | PartialDt
> \PartialDx /
>
> /PartialD
> ApplyFunction(x,t)=uminus0vApplyFunction(x,t) |--------- uApplyFunction(x,t)
> \PartialDx
>
> \ 2 PartialD
> |+vApplyFunction(x,t) uApplyFunction(x,t) ,--------- vApplyFunction(x,t)=
> / PartialDx
>
> ]
> ]
> uminus0uApplyFunction(x,t) vApplyFunction(x,t)]
> ]
> ]
[ 2 ]
sys4 := [u[t] + 2 u u[x] - u[x, x] = 0, v[t] = -v u[x] + v u , v[x] = -u v]
Apart from the solution to the general or main case, nonlinear systems usually have singular solutions, which solve the subsystems appearing during the splitting into cases of the differential elimination process. Depending on the case, the computation of the singular cases may be quite time consuming. To compute only the solution to the general case, one can give the optional argument singsol = false.
> sol:=pdsolveApplyFunction(sys4,[u,v],singsol=false)
/ / 2 \
| (1/2) |/ / (1/2) \\ |
| _c[1] \\exp\_c[1] x// _C1 - _C2/
sol := < u = - -------------------------------------------,
| 2
| / / (1/2) \\
\ \exp\_c[1] x// _C1 + _C2

\
|

/ (1/2) \ _C3 exp(_c[1] t) _C2|
v = exp\_c[1] x/ _C3 exp(_c[1] t) _C1 + -------------------- >
/ (1/2) \ |
exp\_c[1] x/ |
/
> pdetestApplyFunction(sol,sys4)
[0, 0, 0]
Without giving singsol = false, pdsolve attempts to compute the solutions to both the general and singular cases, in this example leading to
> pdsolveApplyFunction(sys4,[u,v])
/
|
/ _C3 \ |
{ u = - ----- - _C2 tanh(_C1 + _C2 x + _C3 t), v = 0 }, <
\ 2 _C2 / |
|
\

/ 2 \
(1/2) |/ / (1/2) \\ |

_c[1] \\exp\_c[1] x// _C1 - _C2/
u = - -------------------------------------------,
2
/ / (1/2) \\
\exp\_c[1] x// _C1 + _C2

\
|

/ (1/2) \ _C3 exp(_c[1] t) _C2|
v = exp\_c[1] x/ _C3 exp(_c[1] t) _C1 + -------------------- >
/ (1/2) \ |
exp\_c[1] x/ |
/
In the singular case above, v = 0, only Burger's equation remains in the system and the solution computed is a Traveling Wave Solution. The introduction of a potential v in the system leads to the general solution only when v is different from zero.
To see the cases, pass to PDEtools[casesplit] the same arguments passed to pdsolve.
> casesplitApplyFunction(sys4,[u,v])
[ v[x] ]
[u = - ----, v[x, x] = v[t]] &where [v <> 0],
[ v ]

[u[x, x] = u[t] + 2 u u[x], v = 0] &where []

Example 5.
Compute the solutions of the following overdetermined linear system of 38 PDEs that arises when computing point symmetries for the wave equation (more like an exercise concerning the solving of a system with many equations). The system is given by
> declareApplyFunction((?,?)ApplyFunction(x,y,z,t,u))
xi(x, y, z, t, u) will now be displayed as xi
eta(x, y, z, t, u) will now be displayed as eta
> [
> [PartialD PartialD
> sys5:=[--------- ?[1]ApplyFunction(x,y,z,t,u)=0,--------- ?[1]ApplyFunction(x,
> [PartialDu PartialDx
> [
>
> /PartialD \ PartialD
> y,z,t,u) - |--------- ?[2]ApplyFunction(x,y,z,t,u)|=0,--------- ?[2]
> \PartialDy / PartialDu
>
> /PartialD \
> ApplyFunction(x,y,z,t,u)=0,uminus0|--------- ?[1]ApplyFunction(x,y,z,t,u)|
> \PartialDy /
>
> /PartialD \ PartialD
> - |--------- ?[2]ApplyFunction(x,y,z,t,u)|=0,--------- ?[3]ApplyFunction(x,
> \PartialDx / PartialDu
>
> PartialD /PartialD
> y,z,t,u)=0,--------- ?[1]ApplyFunction(x,y,z,t,u) - |--------- ?[3]
> PartialDx \PartialDz
>
> \ /PartialD \
> ApplyFunction(x,y,z,t,u)|=0,uminus0|--------- ?[3]ApplyFunction(x,y,z,t,u)|
> / \PartialDy /
>
> /PartialD \ /PartialD
> - |--------- ?[2]ApplyFunction(x,y,z,t,u)|=0,uminus0|--------- ?[1]
> \PartialDz / \PartialDz
>
> \ /PartialD \
> ApplyFunction(x,y,z,t,u)| - |--------- ?[3]ApplyFunction(x,y,z,t,u)|=0,
> / \PartialDx /
>
> PartialD PartialD
> --------- ?[4]ApplyFunction(x,y,z,t,u)=0,--------- ?[3]ApplyFunction(x,y,z,t
> PartialDu PartialDt
>
> /PartialD \ PartialD
> ,u) - |--------- ?[4]ApplyFunction(x,y,z,t,u)|=0,--------- ?[2]ApplyFunction
> \PartialDz / PartialDt
>
> /PartialD \ PartialD
> (x,y,z,t,u) - |--------- ?[4]ApplyFunction(x,y,z,t,u)|=0,--------- ?[1]
> \PartialDy / PartialDt
>
> /PartialD \
> ApplyFunction(x,y,z,t,u) - |--------- ?[4]ApplyFunction(x,y,z,t,u)|=0,
> \PartialDx /
>
> /PartialD \ PartialD
> uminus0|--------- ?[1]ApplyFunction(x,y,z,t,u)|+--------- ?[4]ApplyFunction(
> \PartialDx / PartialDt
>
> 2 2
> PartialD PartialD
> x,y,z,t,u)=0,---------- ?[1]ApplyFunction(x,y,z,t,u)+---------- ?[1]
> 2 2
> PartialDy PartialDz
>
> / 2 \
> |PartialD |
> ApplyFunction(x,y,z,t,u) - |---------- ?[1]ApplyFunction(x,y,z,t,u)|+
> | 2 |
> \PartialDt /
>
> 2 2
> PartialD PartialD
> ---------- ?[1]ApplyFunction(x,y,z,t,u)=0,---------- ?[1]ApplyFunction(x,y,z
> 2 2
> PartialDx PartialDu
>
> 2 2
> PartialD PartialD
> ,t,u)=0,------------------- ?[1]ApplyFunction(x,y,z,t,u)+---------- ?[1]
> PartialDx PartialDu 2
> PartialDx
>
> 2
> PartialD
> ApplyFunction(x,y,z,t,u)=0,------------------- ?[1]ApplyFunction(x,y,z,t,u)+
> PartialDy PartialDx
>
> 2 / 2
> PartialD |PartialD
> ------------------- ?[1]ApplyFunction(x,y,z,t,u)=0,uminus0|---------- ?[1]
> PartialDy PartialDu | 2
> \PartialDy
>
> \ 2
> | PartialD
> ApplyFunction(x,y,z,t,u)|+------------------- ?[1]ApplyFunction(x,y,z,t,u)=0
> | PartialDx PartialDu
> /
>
> 2 2
> PartialD PartialD
> ,------------------- ?[1]ApplyFunction(x,y,z,t,u)+------------------- ?[1]
> PartialDz PartialDx PartialDz PartialDu
>
> 2
> PartialD
> ApplyFunction(x,y,z,t,u)=0,------------------- ?[1]ApplyFunction(x,y,z,t,u)=
> PartialDz PartialDy
>
> / 2 \ 2
> |PartialD | PartialD
> 0,uminus0|---------- ?[1]ApplyFunction(x,y,z,t,u)|+------------------- ?[1]
> | 2 | PartialDx PartialDu
> \PartialDz /
>
> / 2
> | PartialD
> ApplyFunction(x,y,z,t,u)=0,uminus0|------------------- ?[1]ApplyFunction(x,y
> \PartialDu PartialDt
>
> \ / 2 \
> | | PartialD |
> ,z,t,u)| - |------------------- ?[1]ApplyFunction(x,y,z,t,u)|=0,
> / \PartialDx PartialDt /
>
> 2 2
> PartialD PartialD
> ------------------- ?[1]ApplyFunction(x,y,z,t,u)=0,------------------- ?[1]
> PartialDy PartialDt PartialDz PartialDt
>
> 2
> PartialD
> ApplyFunction(x,y,z,t,u)=0,---------- ?[1]ApplyFunction(x,y,z,t,u)+
> 2
> PartialDt
>
> 2 / 2
> PartialD |PartialD
> ------------------- ?[1]ApplyFunction(x,y,z,t,u)=0,uminus0|---------- ?[2]
> PartialDx PartialDu | 2
> \PartialDz
>
> \ 2
> | PartialD
> ApplyFunction(x,y,z,t,u)|+------------------- ?[1]ApplyFunction(x,y,z,t,u)=0
> | PartialDy PartialDu
> /
>
> 2 2
> PartialD PartialD
> ,------------------- ?[2]ApplyFunction(x,y,z,t,u)=0,---------- ?[2]
> PartialDz PartialDt 2
> PartialDt
>
> 2
> PartialD
> ApplyFunction(x,y,z,t,u)+------------------- ?[1]ApplyFunction(x,y,z,t,u)=0,
> PartialDy PartialDu
>
> 2 2
> PartialD PartialD
> ---------- ?[3]ApplyFunction(x,y,z,t,u)+------------------- ?[1]
> 2 PartialDz PartialDu
> PartialDt
>
> 3
> PartialD
> ApplyFunction(x,y,z,t,u)=0,-------------------- ?[1]ApplyFunction(x,y,z,t,u)
> 2
> PartialDx PartialDu
>
> 3
> PartialD
> =0,----------------------------- ?[1]ApplyFunction(x,y,z,t,u)=0,
> PartialDy PartialDx PartialDu
>
> 3
> PartialD
> -------------------- ?[1]ApplyFunction(x,y,z,t,u)=0,
> 2
> PartialDy PartialDu
>
> 3
> PartialD
> ----------------------------- ?[1]ApplyFunction(x,y,z,t,u)=0,
> PartialDz PartialDx PartialDu
>
> 3
> PartialD
> ----------------------------- ?[1]ApplyFunction(x,y,z,t,u)=0,
> PartialDz PartialDy PartialDu
>
> 3
> PartialD
> -------------------- ?[1]ApplyFunction(x,y,z,t,u)=0,
> 2
> PartialDz PartialDu
>
> 3
> PartialD
> ----------------------------- ?[1]ApplyFunction(x,y,z,t,u)=0,
> PartialDx PartialDu PartialDt
>
> 3
> PartialD
> ----------------------------- ?[1]ApplyFunction(x,y,z,t,u)=0,
> PartialDy PartialDu PartialDt
>
> 3 ]
> PartialD ]
> ----------------------------- ?[1]ApplyFunction(x,y,z,t,u)=0]
> PartialDz PartialDu PartialDt ]
> ]
sys5 := [xi[1, u] = 0, xi[1, x] - xi[2, y] = 0, xi[2, u] = 0,

-xi[1, y] - xi[2, x] = 0, xi[3, u] = 0, xi[1, x] - xi[3, z] = 0,

-xi[3, y] - xi[2, z] = 0, -xi[1, z] - xi[3, x] = 0, xi[4, u] = 0,

xi[3, t] - xi[4, z] = 0, xi[2, t] - xi[4, y] = 0, xi[1, t] - xi[4, x] = 0,

-xi[1, x] + xi[4, t] = 0,

eta[1, y, y] + eta[1, z, z] - eta[1, t, t] + eta[1, x, x] = 0,

eta[1, u, u] = 0, eta[1, u, x] + xi[1, x, x] = 0,

xi[1, x, y] + eta[1, u, y] = 0, -xi[1, y, y] + eta[1, u, x] = 0,

xi[1, x, z] + eta[1, u, z] = 0, xi[1, y, z] = 0,

-xi[1, z, z] + eta[1, u, x] = 0, -eta[1, t, u] - xi[1, t, x] = 0,

xi[1, t, y] = 0, xi[1, t, z] = 0, xi[1, t, t] + eta[1, u, x] = 0,

-xi[2, z, z] + eta[1, u, y] = 0, xi[2, t, z] = 0,

xi[2, t, t] + eta[1, u, y] = 0, xi[3, t, t] + eta[1, u, z] = 0,

eta[1, u, x, x] = 0, eta[1, u, x, y] = 0, eta[1, u, y, y] = 0,

eta[1, u, x, z] = 0, eta[1, u, y, z] = 0, eta[1, u, z, z] = 0,

eta[1, t, u, x] = 0, eta[1, t, u, y] = 0, eta[1, t, u, z] = 0]

> nopsApplyFunction(sys5)
38
The solution is given by
> sol:=pdsolveApplyFunction(sys5)
/ /
| 1 |
sol := < eta[1] = ----------------------------------------------------- \_C13
| / (1/2) \ / (1/2) \ / (1/2) \
\ exp\_c[1] x/ exp\_c[2] y/ exp\_c[3] z/

/ 2 \ / 2 \ /

| / / (1/2) \\ | |/ / (1/2) \\ | |
\_C10 \exp\_c[3] z// + _C11/ \\exp\_c[2] y// _C8 + _C9/ \

2 \
/ / (1/2) \\ | / (1/2) \

\exp\_c[1] x// _C6 + _C7/ cos\(-_c[1] - _c[2] - _c[3]) t/ + _C12

/ 2 \ / 2 \ /

| / / (1/2) \\ | |/ / (1/2) \\ | |
\_C10 \exp\_c[3] z// + _C11/ \\exp\_c[2] y// _C8 + _C9/ \

2 \
/ / (1/2) \\ | / (1/2) \ /

\exp\_c[1] x// _C6 + _C7/ sin\(-_c[1] - _c[2] - _c[3]) t/ + u exp\

(1/2) \ / (1/2) \ / (1/2) \

_c[1] x/ exp\_c[2] y/ exp\_c[3] z/ (_C1 t + _C2 x + _C5 y

\
| 1 2

+ _C3 z + _C4)/, xi[1] = - - _C2 x
2

1 1 / 2 2 2\

+ - (-2 _C5 y - 2 _C3 z - 2 _C1 t + 2 _C16) x + - \y + z - t / _C2
2 2

1 2
+ z _C15 + y _C14 + _C17 + t _C18, xi[2] = - - _C5 y

2

1 1 / 2 2 2\

+ - (-2 _C2 x - 2 _C3 z - 2 _C1 t + 2 _C16) y + - \x - t + z / _C5
2 2

1 2
+ _C21 z + _C19 t - x _C14 + _C20, xi[3] = - - _C3 z

2

1 1 / 2 2 2\

+ - (-2 _C5 y + 2 _C16 - 2 _C2 x - 2 _C1 t) z + - \x + y - t / _C3
2 2

1 2
- y _C21 - x _C15 + _C22 t + _C23, xi[4] = - - _C1 t

2

1 1 / 2 2 2\

+ - (-2 _C2 x - 2 _C5 y + 2 _C16 - 2 _C3 z) t + - \-x - y - z / _C1
2 2

\
|

+ z _C22 + y _C19 + x _C18 + _C24 >
|
/
> pdetestApplyFunction(sol,sys5)
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Example 6.
The following is a nonlinear PDE system example in five unknowns, containing four PDEs, extracted from Landau and Lifshitz's "Fluid Mechanics" book. To obtain enhanced display and simplify the input
> PDEtools[declare]ApplyFunction((u,w,?)ApplyFunction(x,z,t),pApplyFunction(x,z)
>
> ,TsApplyFunction(z))
u(x, z, t) will now be displayed as u
w(x, z, t) will now be displayed as w
Tau(x, z, t) will now be displayed as Tau
p(x, z) will now be displayed as p
Ts(z) will now be displayed as Ts
> withApplyFunction(DEtools,diff_table)
[diff_table]
> U:=diff_tableApplyFunction(uApplyFunction(x,z,t));

W:=diff_tableApplyFunction(wApplyFunction(x,z,t));

P:=diff_tableApplyFunction(pApplyFunction(x,z));

T:=diff_tableApplyFunction(?ApplyFunction(x,z,t))

The PDE system is:
> e1:=U[x]+W[z]=0
e1 := u[x] + w[z] = 0
> e2:=U[t]+U[[]] U[x]+W[[]] U[z]=uminus0P[x]+? (U[x,x]+U[z,z])
e2 := u[t] + u u[x] + w u[z] = -p[x] + sigma (u[x, x] + u[z, z])
> e3:=W[t]+U[[]] W[x]+W[[]] W[z]=uminus0P[z]+? (W[x,x]+W[z,z])+? ? (T[[]] - Ts
>
> ApplyFunction(z))
e3 := w[t] + u w[x] + w w[z] = -p[z] + sigma (w[x, x] + w[z, z])

+ sigma rho (Tau - Ts)

> e4:=T[t]+U[[]] T[x]+W[[]] T[z]=T[x,x]+T[z,z]
e4 := Tau[t] + u Tau[x] + w Tau[z] = Tau[x, x] + Tau[z, z]
> sys:=[e1,e2,e3,e4]
This system represents a 2-D (coordinates x, z) fluid layer, open to air, subjected to a vertical temperature gradient (Nabla of T(x,z,t)) and embedded in a gravitational field. The constants sigma and rho are related to the viscosity, thermal diffusivity, the gravitational constant, and the coefficient of thermal expansion. A set of traveling wave solutions can be computed for this system in 2 seconds (typical PC 2 GigaHertz, 2004). The nonsingular case after removing redundant constant solutions:
> pdsolveApplyFunction(sys,[u,w,?,p,Ts],HINT=TWS,singsol=false)
/
| 1 /
{Tau = _C8, Ts = _C8, p = _C9, u = _C5, w = _C7}, < Tau = ----------------- \
| _C2 _C3 sigma rho
\
2 3 2
-_C6 _C2 _C7 _C3 + _C2 _C3 sigma rho _C10 - _C6 _C3 _C7 - _C6 _C2 _C4

3 2 2 \

- _C6 _C2 _C5 - _C6 _C3 _C4 - _C6 _C3 _C5 _C2/, Ts = _C10

/ 4 2 2 4\

2 _C6 \_C2 + 2 _C2 _C3 + _C3 / tanh(_C1 + _C2 x + _C3 z + _C4 t)
+ ------------------------------------------------------------------- -
_C2 _C3 rho

1 // 2 3 2 3 2

----------------- \\_C2 _C7 _C3 + _C3 _C7 + _C2 _C4 + _C2 _C5 + _C3 _C4
_C2 _C3 sigma rho

2 \ 2\
+ _C3 _C5 _C2/ _C6 tanh(_C1 + _C2 x + _C3 z + _C4 t) /

/ 4 2 2 4\ 3
2 _C6 \_C2 + 2 _C2 _C3 + _C3 / tanh(_C1 + _C2 x + _C3 z + _C4 t)
- --------------------------------------------------------------------, p =
_C2 _C3 rho

(_C4 + _C7 _C3 + _C5 _C2) _C6 tanh(_C1 + _C2 x + _C3 z + _C4 t)

_C9 - ---------------------------------------------------------------
_C2

/ 2 2\ 2

\_C2 + _C3 / _C6 sigma tanh(_C1 + _C2 x + _C3 z + _C4 t)
- ----------------------------------------------------------,
_C2

u = _C5 + _C6 tanh(_C1 + _C2 x + _C3 z + _C4 t),

\
_C6 _C2 tanh(_C1 + _C2 x + _C3 z + _C4 t)|
w = _C7 - ----------------------------------------- >
_C3 |
/
For p = constant, this system also admits solutions in terms of the WeierstrassP elliptic function; this is last in the sequence of nonsingular solutions (operand [-1]).
> pdsolveApplyFunction(sys,[u,w,?,p,Ts],HINT=TWSApplyFunction(WeierstrassP),
>
> singsol=false,remove_redundant=false)[uminus01]
/
{ Tau = _C14 + _C11 WeierstrassP(_C3 + I _C5 x + _C5 z + _C6 t, _C2, _C1)
\

2
+ _C12 WeierstrassP(_C3 + I _C5 x + _C5 z + _C6 t, _C2, _C1) , Ts = _C14

+ _C11 WeierstrassP(_C3 + I _C5 x + _C5 z + _C6 t, _C2, _C1)

2
+ _C12 WeierstrassP(_C3 + I _C5 x + _C5 z + _C6 t, _C2, _C1) , p = _C13,

u = _C7 + _C8 WeierstrassP(_C3 + I _C5 x + _C5 z + _C6 t, _C2, _C1), w =
I _C7 _C5 + _C6
- ---------------
_C5

\
- I _C8 WeierstrassP(_C3 + I _C5 x + _C5 z + _C6 t, _C2, _C1) }

/
See Also
casesplit, dchange, declare, DEtools, diff, diffalg, dpolyform, dsolve, dsolve,system, gensys, int, odetest, pdetest, PDEtools, pdsolve, Rif, TWSolutions, solve
References
Cheb-Terrab, E.S. "A Computational Approach for the Exact Solving of Systems of Partial Differential Equations". Submitted to Computer Physics Communications, 2001.