pdsolve
find exact solutions for systems of partial differential equations (PDEs)
Calling Sequence
Parameters
Description
Examples
References
pdsolve(PDE_sys, optional_1, optional_2,...)
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
generalsolution
return only a general solution to the system, or NULL otherwise
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 DifferentialAlgebra package
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. pdsolve also works with anticommutative variables set using the Physics package using the approach explained in PerformOnAnticommutativeSystem.
The solving process
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.
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 g,h 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 DifferentialAlgebra 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[Rif] package. This default can be changed by using the optional argument diffalg, anywhere in the calling sequence, which causes the computation to use the DifferentialAlgebra package.
For further mathematical conventions related to the output by pdsolve, see the conventions used in solve, dsolve, and int.
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:
with(PDEtools, casesplit, declare);
casesplit,declare
with(DEtools,gensys);
gensys
declare( y(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 '
ode[11] := diff(y(x),x,x)+a*x^r*y(x)^n = 0;
ode11≔y''+a⁢xr⁢yn=0
The PDE system satisfied by the symmetries, that is, infinitesimals ξ,η of the symmetry generator, of the ODE above is given by
declare( (xi,eta)(x,y) );
ξ⁡x,y⁢will now be displayed as⁢ξ
η⁡x,y⁢will now be displayed as⁢η
sys := [ gensys( ode[11], [xi,eta](x,y) ) ]:
for _eq in sys do _eq=0 end do;
ξy,y=0
ηy,y−2⁢ξx,y=0
3⁢xr⁢yn⁢ξy⁢a+2⁢ηx,y−ξx,x=0
2⁢ξx⁢xr⁢yn⁢a−xr⁢yn⁢ηy⁢a+η⁢a⁢xr⁢yn⁢ny+ξ⁢a⁢xr⁢r⁢ynx+ηx,x=0
This is a second order linear PDE system, with two unknowns η⁡x,y,ξ⁡x,y and four equations. Its general solution is given by
sol := pdsolve(sys);
sol≔η=−y⁢c__1⁢r+2n−1,ξ=c__1⁢x
Since symmetries are defined up to a multiplicative constant, we can drop c__1. 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).
pdetest(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⁢xr⁢yn=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 := [op(sys), n<>1]:
Next we indicate to pdsolve that n and r are parameters of the problem, thus arriving at the desired solutions.
sol1 := pdsolve(sys1, parameters={n,r});
sol1≔n=2,r=−5,η=y⁢c__1⁢x+3⁢c__2,ξ=x⁢c__1⁢x+c__2,n=2,r=−207,η=−2⁢−6⁢c__1⁢x2−98⁢x87⁢c__1⁢a⁢y−147⁢c__2⁢a⁢x⁢y343⁢x⁢a,ξ=c__1⁢x87+c__2⁢x,n=2,r=−157,η=3⁢x67⁢c__1⁢a⁢y+x⁢c__2⁢a⁢y−12⁢c__1497⁢a⁢x,ξ=c__1⁢x67+c__2⁢x,n=2,r=r,η=−y⁢c__1⁢r+2,ξ=c__1⁢x,n=−r−3,r=r,η=y⁢c__1⁢x+c__2⁢r+4⁢c__1⁢x+2⁢c__2r+4,ξ=x⁢c__1⁢x+c__2,n=n,r=r,η=−y⁢c__1⁢r+2n−1,ξ=c__1⁢x
map(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.
declare( F(r,s), H(r), G(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
sys2 := [-diff(F(r,s),r,r) + diff(F(r,s),s,s) + diff(H(r),r) + diff(G(s),s) + s = 0, diff(F(r,s),r,r) + 2*diff(F(r,s),r,s) + diff(F(r,s),s,s) - diff(H(r),r) + diff(G(s), s)-r = 0];
sys2≔−Fr,r+Fs,s+Hr+Gs+s=0,Fr,r+2⁢Fr,s+Fs,s−Hr+Gs−r=0
The solution for the unknowns G, H is given by the following expression, where determining whether this solution is or is not a general solution, is non-obvious. The userinfo is helpful here to understand the solution process:
infolevel[pdsolve] := 3:
sol := pdsolve(sys2);
-> Solving ordering for the dependent variables of the PDE system: [F(r,s), H(r), G(s)]
-> Solving ordering for the independent variables (can be changed using the ivars option): [r, s]
tackling triangularized subsystem with respect to F(r,s)
First set of solution methods (general or quasi general solution)
-> trying differential factorization for linear PDEs ...
<- differential factorization successful.
<- First set of solution methods successful
tackling triangularized subsystem with respect to H(r)
tackling triangularized subsystem with respect to G(s)
<- Returning a *general* solution
sol≔F=f__1⁡s+f__2⁡r+f__3⁡s−r−r2⁢r−3⁢s12,G=−s24−f__1s+c__2,H=−r24+f__2r+c__1
and we see that the solution given is general, depending on 3 arbitrary functions, f__1⁡s,f__2⁡r,f__3⁡s−r.
infolevel[pdsolve] := 1:
pdetest(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.
sys3 := [-y*diff(f(x,y,z,t),x) + z^2*diff(f(x,y,z,t),z) + 3*t*z*diff(f(x,y,z,t),t) - 3*t^2-4*f(x,y,z,t)*z = 0, -y*diff(f(x,y,z,t),y) - z*diff(f(x,y,z,t),z) - t*diff(f(x,y,z,t),t) + f(x,y,z,t) = 0, -x*diff(f(x,y,z,t),y) - diff(f(x,y,z,t),z) = 0]:
for _eq in sys3 do _eq; end do;
−y⁢fx+z2⁢fz+3⁢t⁢z⁢ft−3⁢t2−4⁢f⁡x,y,z,t⁢z=0
−y⁢fy−z⁢fz−t⁢ft+f⁡x,y,z,t=0
−x⁢fy−fz=0
The solving technique makes use of differential invariants and the solution is given by
sol := pdsolve(sys3);
sol≔f⁡x,y,z,t=3⁢t2⁢xx⁢z−y+t32⁢c__1x⁢z−y
pdetest(sol, sys3);
0,0,0
Example 4.
Compute the solution of the following nonlinear system, consisting of Burgers' equation and a possible potential.
declare(u(x,t), v(x,t));
u⁡x,t⁢will now be displayed as⁢u
v⁡x,t⁢will now be displayed as⁢v
sys4 := [diff(u(x,t),t)+2*u(x,t)*diff(u(x,t),x)-diff(u(x,t),x,x) = 0, diff(v(x,t),t) = -v(x,t)*diff(u(x,t),x)+v(x,t)*u(x,t)^2, diff(v(x,t),x) = -u(x,t)*v(x,t)];
sys4≔ut+2⁢u⁢ux−ux,x=0,vt=−v⁢ux+v⁢u2,vx=−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 := pdsolve(sys4, [u,v], singsol=false);
sol≔u=−_c1⁢ⅇ_c1⁢x2⁢c__1−c__2ⅇ_c1⁢x2⁢c__1+c__2,v=c__3⁢ⅇ_c1⁢t⁢c__1⁢ⅇ_c1⁢x+c__3⁢ⅇ_c1⁢t⁢c__2ⅇ_c1⁢x
pdetest(sol,sys4);
Without giving singsol = false, pdsolve attempts to compute the solutions to both the general and singular cases, in this example leading to
pdsolve(sys4, [u,v]);
* Using tau = tanh(t*C[2]+x*C[1]+C[0])
* Equivalent ODE system: {C[1]^2*(tau^2-1)^2*diff(diff(u(tau),tau),tau)+(2*C[1]^2*(tau^2-1)*tau+2*u(tau)*C[1]*(tau^2-1)+C[2]*(tau^2-1))*diff(u(tau),tau)}
* Power series solution [1]: {u(tau) = tau*A[1,1]+A[1,0]}
u=−c__2⁢tanh⁡c__2⁢x+c__3⁢t+c__1−c__32⁢c__2,v=0,u=−_c1⁢ⅇ_c1⁢x2⁢c__1−c__2ⅇ_c1⁢x2⁢c__1+c__2,v=c__3⁢ⅇ_c1⁢t⁢c__1⁢ⅇ_c1⁢x+c__3⁢ⅇ_c1⁢t⁢c__2ⅇ_c1⁢x
In the singular case above, v=0, only Burgers' 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.
casesplit(sys4,[u,v]);
u=−vxv,vx,x=vtwherev≠0,ux,x=2⁢u⁢ux+ut,v=0where
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
declare( (xi,eta)(x,y,z,t,u) );
ξ⁡x,y,z,t,u⁢will now be displayed as⁢ξ
η⁡x,y,z,t,u⁢will now be displayed as⁢η
sys5 := [diff(xi[1](x,y,z,t,u),u) = 0, diff(xi[1](x,y,z,t,u),x)-diff(xi[2](x,y,z,t,u),y) = 0, diff(xi[2](x,y,z,t,u),u) = 0, -diff(xi[1](x,y,z,t,u),y)-diff(xi[2](x,y,z,t,u),x) = 0, diff(xi[3](x,y,z,t,u),u) = 0, diff(xi[1](x,y,z,t,u),x)-diff(xi[3](x,y,z,t,u),z) = 0, -diff(xi[3](x,y,z,t,u),y)-diff(xi[2](x,y,z,t,u),z) = 0, -diff(xi[1](x,y,z,t,u),z)-diff(xi[3](x,y,z,t,u),x) = 0, diff(xi[4](x,y,z,t,u),u) = 0, diff(xi[3](x,y,z,t,u),t)-diff(xi[4](x,y,z,t,u),z) = 0, diff(xi[2](x,y,z,t,u),t)-diff(xi[4](x,y,z,t,u),y) = 0, diff(xi[1](x,y,z,t,u),t)-diff(xi[4](x,y,z,t,u),x) = 0, -diff(xi[1](x,y,z,t,u),x)+diff(xi[4](x,y,z,t,u),t) = 0, diff(eta[1](x,y,z,t,u),y,y)+diff(eta[1](x,y,z,t,u),z,z)-diff(eta[1](x,y,z,t, u),t,t)+diff(eta[1](x,y,z,t,u),x,x) = 0, diff(eta[1](x,y,z,t,u),u,u) = 0, diff(eta[1](x,y,z,t,u),u,x)+diff(xi[1](x,y,z,t,u),x,x) = 0, diff(xi[1](x,y,z,t,u),x,y)+diff(eta[1](x,y,z,t,u),u,y) = 0, -diff(xi[1](x,y,z,t,u),y,y)+diff(eta[1](x,y,z,t,u),u,x) = 0, diff(xi[1](x,y,z,t,u),x,z)+diff(eta[1](x,y,z,t,u),u,z) = 0, diff(xi[1](x,y,z,t,u),y,z) = 0, -diff(xi[1](x,y,z,t,u),z,z)+diff(eta[1](x,y,z,t,u),u,x) = 0, -diff(eta[1](x,y,z,t,u),t,u)-diff(xi[1](x,y,z,t,u),t,x) = 0, diff(xi[1](x,y,z,t,u),t,y) = 0, diff(xi[1](x,y,z,t,u),t,z) = 0, diff(xi[1](x,y,z,t,u),t,t)+diff(eta[1](x,y,z,t,u),u,x) = 0, -diff(xi[2](x,y,z,t,u),z,z)+diff(eta[1](x,y,z,t,u),u,y) = 0, diff(xi[2](x,y,z,t,u),t,z) = 0, diff(xi[2](x,y,z,t,u),t,t)+diff(eta[1](x,y,z,t,u),u,y) = 0, diff(xi[3](x,y,z,t,u),t,t)+diff(eta[1](x,y,z,t,u),u,z) = 0, diff(eta[1](x,y,z,t,u),u,x,x) = 0, diff(eta[1](x,y,z,t,u),u,x,y) = 0, diff(eta[1](x,y,z,t,u),u,y,y) = 0, diff(eta[1](x,y,z,t,u),u,x,z) = 0, diff(eta[1](x,y,z,t,u),u,y,z) = 0, diff(eta[1](x,y,z,t,u),u,z,z) = 0, diff(eta[1](x,y,z,t,u),t,u,x) = 0, diff(eta[1](x,y,z,t,u),t,u,y) = 0, diff(eta[1](x,y,z,t,u),t,u,z) = 0];
sys5≔ξ1u=0,ξ1x−ξ2y=0,ξ2u=0,−ξ1y−ξ2x=0,ξ3u=0,ξ1x−ξ3z=0,−ξ3y−ξ2z=0,−ξ1z−ξ3x=0,ξ4u=0,ξ3t−ξ4z=0,ξ2t−ξ4y=0,ξ1t−ξ4x=0,−ξ1x+ξ4t=0,η1y,y+η1z,z−η1t,t+η1x,x=0,η1u,u=0,η1u,x+ξ1x,x=0,ξ1x,y+η1u,y=0,−ξ1y,y+η1u,x=0,ξ1x,z+η1u,z=0,ξ1y,z=0,−ξ1z,z+η1u,x=0,−η1t,u−ξ1t,x=0,ξ1t,y=0,ξ1t,z=0,ξ1t,t+η1u,x=0,−ξ2z,z+η1u,y=0,ξ2t,z=0,ξ2t,t+η1u,y=0,ξ3t,t+η1u,z=0,η1u,x,x=0,η1u,x,y=0,η1u,y,y=0,η1u,x,z=0,η1u,y,z=0,η1u,z,z=0,η1t,u,x=0,η1t,u,y=0,η1t,u,z=0
nops(sys5);
38
The solution is given by
sol := pdsolve(sys5);
sol≔η1=c__13⁢c__10⁢ⅇ_c3⁢z2+c__11⁢c__8⁢ⅇ_c2⁢y2+c__9⁢c__6⁢ⅇ_c1⁢x2+c__7⁢cos⁡−_c1−_c2−_c3⁢t+c__12⁢c__10⁢ⅇ_c3⁢z2+c__11⁢c__8⁢ⅇ_c2⁢y2+c__9⁢c__6⁢ⅇ_c1⁢x2+c__7⁢sin⁡−_c1−_c2−_c3⁢t+u⁢ⅇ_c1⁢x⁢ⅇ_c2⁢y⁢ⅇ_c3⁢z⁢c__1⁢t+c__2⁢x+c__3⁢y+c__4⁢z+c__5ⅇ_c3⁢z⁢ⅇ_c2⁢y⁢ⅇ_c1⁢x,ξ1=−c__2⁢x22+−c__1⁢t−c__3⁢y−c__4⁢z+c__17⁢x+−t2+y2+z2⁢c__22+c__16⁢t+c__15⁢z+c__14⁢y+c__18,ξ2=−c__3⁢y22+−c__1⁢t−c__2⁢x−c__4⁢z+c__17⁢y+−t2+x2+z2⁢c__32+c__20⁢t+c__19⁢z−c__14⁢x+c__21,ξ3=−c__4⁢z22+−c__1⁢t−c__2⁢x−c__3⁢y+c__17⁢z+−t2+x2+y2⁢c__42+c__22⁢t−c__19⁢y−c__15⁢x+c__23,ξ4=−c__1⁢t22+−c__2⁢x−c__3⁢y−c__4⁢z+c__17⁢t+−x2−y2−z2⁢c__12+c__20⁢y+c__22⁢z+c__16⁢x+c__24
pdetest(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]((u,w,Tau)(x,z,t),p(x,z),Ts(z));
u⁡x,z,t⁢will now be displayed as⁢u
w⁡x,z,t⁢will now be displayed as⁢w
Τ⁡x,z,t⁢will now be displayed as⁢Τ
p⁡x,z⁢will now be displayed as⁢p
Ts⁡z⁢will now be displayed as⁢Ts
with(DEtools,diff_table);
diff_table
U := diff_table(u(x,z,t)): W := diff_table(w(x,z,t)):
P := diff_table(p(x,z)): T := diff_table(Tau(x,z,t)):
The PDE system is:
e1 := U[x]+W[z] = 0;
e1≔ux+wz=0
e2 := U[t] + U[]*U[x] + W[]*U[z] = -P[x] + sigma*(U[x,x]+U[z,z]);
e2≔u⁢ux+uz⁢w+ut=−px+σ⁢ux,x+uz,z
e3 := W[t]+U[]*W[x]+W[]*W[z] = -P[z]+sigma*(W[x,x]+W[z,z])+sigma*rho*(T[]-Ts(z));
e3≔u⁢wx+w⁢wz+wt=−pz+σ⁢wx,x+wz,z+σ⁢ρ⁢Τ−Ts
e4 := T[t]+U[]*T[x]+W[]*T[z] = T[x,x]+T[z,z];
e4≔Τx⁢u+Τz⁢w+Τt=Τx,x+Τ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:
[pdsolve(sys, [u,w,Tau,p,Ts], HINT = TWS, singsol = false)];
* Using tau = tanh(t*C[3]+x*C[1]+z*C[2]+C[0])
* Equivalent ODE system: {(-C[1]*(tau^2-1)*u(tau)-C[2]*(tau^2-1)*w(tau)-C[3]*(tau^2-1)-2*C[1]^2*(tau^2-1)*tau-2*C[2]^2*(tau^2-1)*tau)*diff(Tau(tau),tau)+(-C[1]^2*(tau^2-1)^2-C[2]^2*(tau^2-1)^2)*diff(diff(Tau(tau),tau),tau), -C[1]*(tau^2-1)*diff(u(tau),tau)-C[2]*(tau^2-1)*diff(w(tau),tau), -sigma*(C[1]^2*(tau^2-1)^2+C[2]^2*(tau^2-1)^2)*diff(diff(u(tau),tau),tau)-C[1]*(tau^2-1)*diff(p(tau),tau)+(-C[1]*(tau^2-1)*u(tau)-C[2]*(tau^2-1)*w(tau)-C[3]*(tau^2-1)-sigma*(2*C[1]^2*(tau^2-1)*tau+2*C[2]^2*(tau^2-1)*tau))*diff(u(tau),tau), -sigma*(C[1]^2*(tau^2-1)^2+C[2]^2*(tau^2-1)^2)*diff(diff(w(tau),tau),tau)-C[2]*(tau^2-1)*diff(p(tau),tau)+(-C[1]*(tau^2-1)*u(tau)-C[2]*(tau^2-1)*w(tau)-C[3]*(tau^2-1)-sigma*(2*C[1]^2*(tau^2-1)*tau+2*C[2]^2*(tau^2-1)*tau))*diff(w(tau),tau)-sigma*rho*(Tau(tau)-Ts(tau))}
* Power series solution [1]: {Tau(tau) = tau*A[3,1]+A[3,0], Ts(tau) = tau^3*A[5,3]+tau^2*A[5,2]+tau*A[5,1]+A[5,0], p(tau) = tau^2*A[4,2]+tau*A[4,1]+A[4,0], u(tau) = tau*A[1,1]+A[1,0], w(tau) = tau*A[2,1]+A[2,0]}
* Power series solution [2]: {Tau(tau) = tau^3*A[3,3]+tau^2*A[3,2]+tau*A[3,1]+A[3,0], Ts(tau) = tau*A[5,1]+A[5,0], p(tau) = tau^2*A[4,2]+tau*A[4,1]+A[4,0], u(tau) = tau*A[1,1]+A[1,0], w(tau) = tau*A[2,1]+A[2,0]}
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]).
pdsolve( sys, [u,w,Tau,p,Ts], HINT = TWS(WeierstrassP), singsol = false, remove_redundant = false)[-1];
* Using tau = WeierstrassP(t*C[3]+x*C[1]+z*C[2]+C[0])
* Equivalent ODE system: {diff(u(tau),tau)*C[1]+diff(w(tau),tau)*C[2], (C[1]*diff(p(tau),tau)+(u(tau)*C[1]+w(tau)*C[2]+C[3])*diff(u(tau),tau))*(4*tau^3-tau*C[-1]-C[-2])^(1/2)-sigma*(4*tau^3-tau*C[-1]-C[-2])*(C[1]^2+C[2]^2)*diff(diff(u(tau),tau),tau)-1/2*sigma*(12*tau^2-C[-1])*(C[1]^2+C[2]^2)*diff(u(tau),tau), (u(tau)*C[1]+w(tau)*C[2]+C[3])*diff(Tau(tau),tau)*(4*tau^3-tau*C[-1]-C[-2])^(1/2)+(-6*tau^2*C[1]^2+1/2*C[-1]*C[1]^2-6*tau^2*C[2]^2+1/2*C[-1]*C[2]^2)*diff(Tau(tau),tau)+(-4*tau^3*C[1]^2-4*tau^3*C[2]^2+tau*C[-1]*C[1]^2+tau*C[-1]*C[2]^2+C[-2]*C[1]^2+C[-2]*C[2]^2)*diff(diff(Tau(tau),tau),tau), (C[2]*diff(p(tau),tau)+(u(tau)*C[1]+w(tau)*C[2]+C[3])*diff(w(tau),tau))*(4*tau^3-tau*C[-1]-C[-2])^(1/2)-sigma*(4*tau^3-tau*C[-1]-C[-2])*(C[1]^2+C[2]^2)*diff(diff(w(tau),tau),tau)-1/2*sigma*(12*tau^2-C[-1])*(C[1]^2+C[2]^2)*diff(w(tau),tau)-sigma*rho*(Tau(tau)-Ts(tau))}
* Power series solution [2]: {Tau(tau) = tau^2*A[3,2]+tau*A[3,1]+A[3,0], Ts(tau) = tau^2*A[5,2]+tau*A[5,1]+A[5,0], p(tau) = tau*A[4,1]+A[4,0], u(tau) = tau*A[1,1]+A[1,0], w(tau) = tau*A[2,1]+A[2,0]}
* Power series solution [3]: {Tau(tau) = tau^3*A[3,3]+tau^2*A[3,2]+tau*A[3,1]+A[3,0], Ts(tau) = tau*A[5,1]+A[5,0], p(tau) = tau^2*A[4,2]+tau*A[4,1]+A[4,0], u(tau) = tau*A[1,1]+A[1,0], w(tau) = tau*A[2,1]+A[2,0]}
* Power series solution [4]: {Tau(tau) = tau*A[3,1]+A[3,0], Ts(tau) = tau^3*A[5,3]+tau^2*A[5,2]+tau*A[5,1]+A[5,0], p(tau) = tau*A[4,1]+A[4,0], u(tau) = tau*A[1,1]+A[1,0], w(tau) = tau*A[2,1]+A[2,0]}
* Power series solution [5]: {Tau(tau) = tau*A[3,1]+A[3,0], Ts(tau) = tau*A[5,1]+A[5,0], p(tau) = tau^2*A[4,2]+tau*A[4,1]+A[4,0], u(tau) = tau*A[1,1]+A[1,0], w(tau) = tau*A[2,1]+A[2,0]}
* Power series solution [6]: {Tau(tau) = tau*A[3,1]+A[3,0], Ts(tau) = tau^2*A[5,2]+tau*A[5,1]+A[5,0], p(tau) = tau^2*A[4,2]+tau*A[4,1]+A[4,0], u(tau) = tau*A[1,1]+A[1,0], w(tau) = tau*A[2,1]+A[2,0]}
* Power series solution [7]: {Tau(tau) = tau^2*A[3,2]+tau*A[3,1]+A[3,0], Ts(tau) = tau*A[5,1]+A[5,0], p(tau) = tau^2*A[4,2]+tau*A[4,1]+A[4,0], u(tau) = tau*A[1,1]+A[1,0], w(tau) = tau*A[2,1]+A[2,0]}
* Power series solution [8]: {Tau(tau) = tau^3*A[3,3]+tau^2*A[3,2]+tau*A[3,1]+A[3,0], Ts(tau) = tau^2*A[5,2]+tau*A[5,1]+A[5,0], p(tau) = tau*A[4,1]+A[4,0], u(tau) = tau*A[1,1]+A[1,0], w(tau) = tau*A[2,1]+A[2,0]}
Example 7.
Both dsolve and pdsolve can now solve PDEs that involve anticommutative variables set using the Physics package using the approach explained in PerformOnAnticommutativeSystem.
with(Physics, Setup);
Setup
Set first θ and Q as suffixes for variables of type/anticommutative (see Setup)
Setup(anticommutativepre = {theta, Q});
* Partial match of 'anticommutativepre' against keyword 'anticommutativeprefix'
_______________________________________________________
anticommutativeprefix=Q,θ
A PDE system example with two unknown anticommutative functions of four variables, two commutative and two anticommutative; to avoid redundant typing in the input that follows and redundant display of information on the screen let's use PDEtools:-diff_table and PDEtools:-declare
PDEtools:-declare(Q(x, y, theta[1], theta[2]));
Q⁡x,y,θ1,θ2⁢will now be displayed as⁢Q
q := PDEtools:-diff_table(Q(x, y, theta[1], theta[2])):
Now we can enter derivatives directly as the function's name indexed by the differentiation variables and see the display the same way; two PDEs
pde[1] := q[x, y, theta[1]] + q[x, y, theta[2]] - q[y, theta[1], theta[2]] = 0;
pde1≔Qx,y,θ1+Qx,y,θ2−Qy,θ1,θ2=0
pde[2] := q[theta[1]] = 0;
pde2≔Qθ1=0
The solution to this system:
pdsolve([pde[1], pde[2]]);
Q=f__1⁡x,y⁢_λ1+f__6⁡x+f__5⁡y⁢θ2
This solution involves an anticommutative constant _lambda2, analogous to the commutative constants _Cn where n is an integer. The arbitrary functions introduced are all commutative, as usual and the Grasmannian parity of the left and right-hand sides is preserved:
Physics:-GrassmannParity((42));
1=1
Cheb-Terrab, E.S. "A Computational Approach for the Exact Solving of Systems of Partial Differential Equations". Submitted to Computer Physics Communications, 2001.
See Also
casesplit
dchange
declare
DEtools
diff
DifferentialAlgebra
dpolyform
dsolve
dsolve,system
int
odetest
pdetest
PDEtools
PerformOnAnticommutativeSystem
Physics
Rif
solve
TWSolutions
Download Help Document