Items tagged with math math Tagged Items Feed

Hi,
The latest update to the differential equations Maple libraries (this week, can be downloaded from the Maplesoft R&D webpage for Differential Equations and Mathematical functions) includes new functionality in pdsolve, regarding whether the solution for a PDE or PDE system is or not a general solution.

In brief, a general solution of a PDE in 1 unknown, that has differential order N, and where the unknown depends on M independent variables, involves N arbitrary functions of M-1 arguments. It is not entirely evident how to extend this definition in the case of a coupled, possibly nonlinear PDE system. However, using differential algebra techniques (automatically used by pdsolve when tackling a PDE system), that extension to define a general solution for a DE system is possible, and also when the system involves ODEs and PDEs, and/or algebraic (that is, non-differential) equations, and/or inequations of the form algebraic*expression <> 0 involving the unknowns, and all of this in the presence of mathematical functions (based on the use of Maple's PDEtools:-dpolyform). This is a very nice case were many different advanced developments come together to naturally solve a problem that otherwise would be rather difficult.

The issues at the center of this Maple development/post are then:

        a) How do you know whether a PDE or PDE system solution returned is a general solution?

        b) How could you indicate to pdsolve that you are only interested in a general PDE or PDE system solution?

The answer to a) is now always present in the last line of the userinfo. So input infolevel[pdsolve] := 3 before calling pdsolve, and check what the last line of the userinfo displayed tells.


The answer to b) is a new option, generalsolution, implemented in pdsolve so that it either returns a general solution or otherwise it returns NULL. If you do not use this new option, then pdsolve works as always: first it tries to compute a general solution and if it fails in doing that it tries to compute a particular solution by separating the variables in different ways, or computing a traveling wave solution or etc. (a number of other well known methods).

 

The examples that follow are from the help page pdsolve,system, and show both the new userinfo telling whether the solution returned is a general one and the option generalsolution at work.The examples are all of differential equation systems but the same userinfos and generalsolution option work as well in the case of a single PDE.

 

 

Example 1.

Solve the determining PDE system for the infinitesimals of the symmetry generator of example 11 from Kamke's book . Tell whether the solution computed is or not a general solution.

infolevel[pdsolve] := 3

3

(1.1)

The PDE system satisfied by the symmetries of Kamke's ODE example number 11 is

sys__1 := [diff(xi(x, y), y, y) = 0, diff(eta(x, y), y, y)-2*(diff(xi(x, y), y, x)) = 0, 3*x^r*y^n*(diff(xi(x, y), y))*a+2*(diff(eta(x, y), y, x))-(diff(xi(x, y), x, x)) = 0, 2*(diff(xi(x, y), x))*x^r*y^n*a-x^r*y^n*(diff(eta(x, y), y))*a+eta(x, y)*a*x^r*y^n*n/y+xi(x, y)*a*x^r*r*y^n/x+diff(eta(x, y), x, x) = 0]

This is a second order linear PDE system, with two unknowns {eta(x, y), xi(x, y)} and four equations. Its general solution is given by the following, where we now can tell that the solution is a general one by reading the last line of the userinfo. Note that because the system is overdetermined, a general solution in this case does not involve any arbitrary function

sol__1 := pdsolve(sys__1)

-> Solving ordering for the dependent variables of the PDE system: [xi(x,y), eta(x,y)]

-> Solving ordering for the independent variables (can be changed using the ivars option): [x, y]
tackling triangularized subsystem with respect to xi(x,y)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
tackling triangularized subsystem with respect to eta(x,y)
<- Returning a *general* solution

 

{eta(x, y) = -_C1*y*(r+2)/(n-1), xi(x, y) = _C1*x}

(1.2)

Next we indicate to pdsolve that n and r are parameters of the problem, and that we want a solution for n <> 1, making more difficult to identify by eye whether the solution returned is or not a general one. Again the last line of the userinfo tells that pdsolve's solution is indeed a general one

`sys__1.1` := [op(sys__1), n <> 1]

[diff(diff(xi(x, y), y), y) = 0, diff(diff(eta(x, y), y), y)-2*(diff(diff(xi(x, y), x), y)) = 0, 3*x^r*y^n*(diff(xi(x, y), y))*a+2*(diff(diff(eta(x, y), x), y))-(diff(diff(xi(x, y), x), x)) = 0, 2*(diff(xi(x, y), x))*x^r*y^n*a-x^r*y^n*(diff(eta(x, y), y))*a+eta(x, y)*a*x^r*y^n*n/y+xi(x, y)*a*x^r*r*y^n/x+diff(diff(eta(x, y), x), x) = 0, n <> 1]

(1.3)

`sol__1.1` := pdsolve(`sys__1.1`, parameters = {n, r})

-> Solving ordering for the dependent variables of the PDE system: [r, n, xi(x,y), eta(x,y)]

-> Solving ordering for the independent variables (can be changed using the ivars option): [x, y]
tackling triangularized subsystem with respect to r
tackling triangularized subsystem with respect to n
tackling triangularized subsystem with respect to xi(x,y)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
tackling triangularized subsystem with respect to eta(x,y)
tackling triangularized subsystem with respect to r
tackling triangularized subsystem with respect to n
tackling triangularized subsystem with respect to xi(x,y)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
tackling triangularized subsystem with respect to eta(x,y)
tackling triangularized subsystem with respect to r
tackling triangularized subsystem with respect to n
tackling triangularized subsystem with respect to xi(x,y)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
tackling triangularized subsystem with respect to eta(x,y)
tackling triangularized subsystem with respect to n
tackling triangularized subsystem with respect to xi(x,y)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
tackling triangularized subsystem with respect to eta(x,y)
tackling triangularized subsystem with respect to n
tackling triangularized subsystem with respect to xi(x,y)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
tackling triangularized subsystem with respect to eta(x,y)
tackling triangularized subsystem with respect to xi(x,y)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
tackling triangularized subsystem with respect to eta(x,y)
<- Returning a *general* solution

 

{n = 2, r = -5, eta(x, y) = y*(_C1*x+3*_C2), xi(x, y) = x*(_C1*x+_C2)}, {n = 2, r = -20/7, eta(x, y) = -(2/343)*(-6*_C1*x^2-98*x^(8/7)*_C1*a*y-147*_C2*a*x*y)/(x*a), xi(x, y) = _C1*x^(8/7)+_C2*x}, {n = 2, r = -15/7, eta(x, y) = -(1/343)*(-49*_C2*a*x*y-147*x^(6/7)*_C1*a*y+12*_C1*x)/(x*a), xi(x, y) = _C1*x^(6/7)+_C2*x}, {n = 2, r = r, eta(x, y) = -_C1*y*(r+2), xi(x, y) = _C1*x}, {n = -r-3, r = r, eta(x, y) = ((_C1*x+_C2)*r+4*_C1*x+2*_C2)*y/(r+4), xi(x, y) = x*(_C1*x+_C2)}, {n = n, r = r, eta(x, y) = -_C1*y*(r+2)/(n-1), xi(x, y) = _C1*x}

(1.4)

map(pdetest, [`sol__1.1`], `sys__1.1`)

[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]

(1.5)

 

Example 2.

Compute the 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.

sys__2 := [-(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]

The solution for the unknowns G, H, is given by the following expression, were again determining whether this solution, that depends on 3 arbitrary functions, _F1(s), _F2(r), _F3(s-r), is or not a general solution, is non-obvious.

sol__2 := pdsolve(sys__2)

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

 

{F(r, s) = _F1(s)+_F2(r)+_F3(s-r)-(1/12)*r^2*(r-3*s), G(s) = -(diff(_F1(s), s))-(1/4)*s^2+_C2, H(r) = diff(_F2(r), r)-(1/4)*r^2+_C1}

(1.6)

pdetest(sol__2, sys__2)

[0, 0]

(1.7)

Example 3.

Compute the solution of the following nonlinear system, consisting of Burger's equation and a possible potential.

sys__3 := [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)]

We see that in this case the solution returned is not a general solution but two particular ones; again the information is in the last line of the userinfo displayed

sol__3 := pdsolve(sys__3, [u, v])

-> Solving ordering for the dependent variables of the PDE system: [v(x,t), u(x,t)]

-> Solving ordering for the independent variables (can be changed using the ivars option): [x, t]
tackling triangularized subsystem with respect to v(x,t)
tackling triangularized subsystem with respect to u(x,t)
First set of solution methods (general or quasi general solution)
Second set of solution methods (complete solutions)
Trying methods for second order PDEs
Third set of solution methods (simple HINTs for separating variables)
PDE linear in highest derivatives - trying a separation of variables by *
HINT = *
Fourth set of solution methods
Trying methods for second order linear PDEs
Preparing a solution HINT ...
Trying HINT = _F1(x)*_F2(t)
Fourth set of solution methods
Preparing a solution HINT ...
Trying HINT = _F1(x)+_F2(t)
Trying travelling wave solutions as power series in tanh ...
* 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)}
* Ordering for functions: [u(tau)]
* Cases for the upper bounds: [[n[1] = 1]]
* Power series solution [1]: {u(tau) = tau*A[1,1]+A[1,0]}
* Solution [1] for {A[i, j], C[k]}: [[A[1,1] = 0], [A[1,0] = -1/2*C[2]/C[1], A[1,1] = -C[1]]]
travelling wave solutions successful.
tackling triangularized subsystem with respect to v(x,t)
First set of solution methods (general or quasi general solution)
Trying differential factorization for linear PDEs ...
Trying methods for PDEs "missing the dependent variable" ...
Second set of solution methods (complete solutions)
Trying methods for second order PDEs
Third set of solution methods (simple HINTs for separating variables)
PDE linear in highest derivatives - trying a separation of variables by *
HINT = *
Fourth set of solution methods
Trying methods for second order linear PDEs
Preparing a solution HINT ...
Trying HINT = _F1(x)*_F2(t)
Third set of solution methods successful
tackling triangularized subsystem with respect to u(x,t)
<- Returning a solution that *is not the most general one*

 

{u(x, t) = -_C2*tanh(_C2*x+_C3*t+_C1)-(1/2)*_C3/_C2, v(x, t) = 0}, {u(x, t) = -_c[1]^(1/2)*((exp(_c[1]^(1/2)*x))^2*_C1-_C2)/((exp(_c[1]^(1/2)*x))^2*_C1+_C2), v(x, t) = _C3*exp(_c[1]*t)*_C1*exp(_c[1]^(1/2)*x)+_C3*exp(_c[1]*t)*_C2/exp(_c[1]^(1/2)*x)}

(1.8)

pdetest(sol__3, sys__3)

[0, 0, 0]

(1.9)

This example is also good for illustrating the other related new feature: one can now request to pdsolve to only compute a general solution (it will return NULL if it cannot achieve that). Turn OFF userinfos and try with this example

infolevel[pdsolve] := 1

This returns NULL:

pdsolve(sys__3, [u, v], generalsolution)

Example 4.

Another where the solution returned is particular, this time for a linear system, conformed by 38 PDEs, also from differential equation symmetry analysis

sys__4 := [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]

There are 38 coupled equations

nops(sys__4)

38

(1.10)

When requesting a general solution pdsolve returns NULL:

pdsolve(sys__4, generalsolution)

A solution that is not a general one, is however computed by default if calling pdsolve without the generalsolution option. In this case again the last line of the userinfo tells that the solution returned is not a general solution

infolevel[pdsolve] := 3

3

(1.11)

sol__4 := pdsolve(sys__4)

-> Solving ordering for the dependent variables of the PDE system: [eta[1](x,y,z,t,u), xi[1](x,y,z,t,u), xi[2](x,y,z,t,u), xi[3](x,y,z,t,u), xi[4](x,y,z,t,u)]

-> Solving ordering for the independent variables (can be changed using the ivars option): [t, x, y, z, u]
tackling triangularized subsystem with respect to eta[1](x,y,z,t,u)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
-> Solving ordering for the dependent variables of the PDE system: [_F1(x,y,z,t), _F2(x,y,z,t)]
-> Solving ordering for the independent variables (can be changed using the ivars option): [t, x, y, z, u]
tackling triangularized subsystem with respect to _F1(x,y,z,t)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
-> Solving ordering for the dependent variables of the PDE system: [_F3(x,y,z), _F4(x,y,z)]
-> Solving ordering for the independent variables (can be changed using the ivars option): [x, y, z, t]
tackling triangularized subsystem with respect to _F3(x,y,z)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
tackling triangularized subsystem with respect to _F4(x,y,z)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
-> Solving ordering for the dependent variables of the PDE system: [_F5(y,z), _F6(y,z)]
-> Solving ordering for the independent variables (can be changed using the ivars option): [y, z, x]
tackling triangularized subsystem with respect to _F5(y,z)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
tackling triangularized subsystem with respect to _F6(y,z)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
-> Solving ordering for the dependent variables of the PDE system: [_F7(z), _F8(z)]
-> Solving ordering for the independent variables (can be changed using the ivars option): [z, y]
tackling triangularized subsystem with respect to _F7(z)
tackling triangularized subsystem with respect to _F8(z)
tackling triangularized subsystem with respect to _F2(x,y,z,t)
First set of solution methods (general or quasi general solution)
Trying differential factorization for linear PDEs ...
Trying methods for PDEs "missing the dependent variable" ...
Second set of solution methods (complete solutions)
Third set of solution methods (simple HINTs for separating variables)
PDE linear in highest derivatives - trying a separation of variables by *
HINT = *
Fourth set of solution methods
Preparing a solution HINT ...
Trying HINT = _F3(x)*_F4(y)*_F5(z)*_F6(t)
Third set of solution methods successful

tackling triangularized subsystem with respect to xi[1](x,y,z,t,u)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
-> Solving ordering for the dependent variables of the PDE system: [_F1(x,z,t), _F2(x,z,t)]
-> Solving ordering for the independent variables (can be changed using the ivars option): [t, x, z, y]
tackling triangularized subsystem with respect to _F1(x,z,t)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
tackling triangularized subsystem with respect to _F2(x,z,t)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful

-> Solving ordering for the dependent variables of the PDE system: [_F3(x,t), _F4(x,t)]
-> Solving ordering for the independent variables (can be changed using the ivars option): [t, x, z]
tackling triangularized subsystem with respect to _F3(x,t)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
tackling triangularized subsystem with respect to _F4(x,t)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
-> Solving ordering for the dependent variables of the PDE system: [_F5(x), _F6(x)]
-> Solving ordering for the independent variables (can be changed using the ivars option): [x, t]
tackling triangularized subsystem with respect to _F5(x)
tackling triangularized subsystem with respect to _F6(x)
tackling triangularized subsystem with respect to xi[2](x,y,z,t,u)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
-> Solving ordering for the dependent variables of the PDE system: [_F1(t), _F2(t)]
-> Solving ordering for the independent variables (can be changed using the ivars option): [t, z]
tackling triangularized subsystem with respect to _F1(t)
tackling triangularized subsystem with respect to _F2(t)
tackling triangularized subsystem with respect to xi[3](x,y,z,t,u)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
tackling triangularized subsystem with respect to xi[4](x,y,z,t,u)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
First set of solution methods successful
<- Returning a solution that *is not the most general one*

 

{eta[1](x, y, z, t, u) = (_C13*(_C10*(exp(_c[3]^(1/2)*z))^2+_C11)*(_C8*(exp(_c[2]^(1/2)*y))^2+_C9)*(_C6*(exp(_c[1]^(1/2)*x))^2+_C7)*cos((-_c[1]-_c[2]-_c[3])^(1/2)*t)+_C12*(_C10*(exp(_c[3]^(1/2)*z))^2+_C11)*(_C8*(exp(_c[2]^(1/2)*y))^2+_C9)*(_C6*(exp(_c[1]^(1/2)*x))^2+_C7)*sin((-_c[1]-_c[2]-_c[3])^(1/2)*t)+u*exp(_c[1]^(1/2)*x)*exp(_c[2]^(1/2)*y)*exp(_c[3]^(1/2)*z)*(_C1*t+_C2*x+_C3*y+_C4*z+_C5))/(exp(_c[1]^(1/2)*x)*exp(_c[2]^(1/2)*y)*exp(_c[3]^(1/2)*z)), xi[1](x, y, z, t, u) = -(1/2)*_C2*x^2+(1/2)*(-2*_C1*t-2*_C3*y-2*_C4*z+2*_C17)*x+(1/2)*(-t^2+y^2+z^2)*_C2+_C16*t+_C15*z+_C14*y+_C18, xi[2](x, y, z, t, u) = -(1/2)*_C3*y^2+(1/2)*(-2*_C1*t-2*_C2*x-2*_C4*z+2*_C17)*y+(1/2)*(-t^2+x^2+z^2)*_C3+_C20*t+_C19*z-_C14*x+_C21, xi[3](x, y, z, t, u) = -(1/2)*_C4*z^2+(1/2)*(-2*_C1*t-2*_C2*x-2*_C3*y+2*_C17)*z+(1/2)*(-t^2+x^2+y^2)*_C4+_C22*t-_C19*y-_C15*x+_C23, xi[4](x, y, z, t, u) = -(1/2)*_C1*t^2+(1/2)*(-2*_C2*x-2*_C3*y-2*_C4*z+2*_C17)*t+(1/2)*(-x^2-y^2-z^2)*_C1+_C20*y+_C22*z+_C16*x+_C24}

(1.12)

pdetest(sol__4, sys__4)

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

(1.13)

Example 5.

Finally, the new userinfos also tell whether a solution is or not a general solution when working with PDEs that involve anticommutative variables  set using the Physics  package

with(Physics, Setup)

[Setup]

(1.14)

Set first theta and Q as suffixes for variables of type/anticommutative  (see Setup )

Setup(anticommutativepre = {Q, theta})

`* Partial match of  'anticommutativepre' against keyword 'anticommutativeprefix'`

 

[anticommutativeprefix = {Q, _lambda, theta}]

(1.15)

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   PDEtools:-declare

PDEtools:-declare(Q(x, y, theta[1], theta[2]))

Q(x, y, theta[1], theta[2])*`will now be displayed as`*Q

(1.16)

q := PDEtools:-diff_table(Q(x, y, theta[1], theta[2]))

table( [(  ) = Q(x, y, theta[1], theta[2]) ] )

(1.17)

Consider the system formed by these two PDEs (because of the q diff_table just defined, we can enter derivatives directly using the function's name indexed by the differentiation variables)

pde[1] := q[x, y, theta[1]]+q[x, y, theta[2]]-q[y, theta[1], theta[2]] = 0

Physics:-diff(diff(diff(Q(x, y, theta[1], theta[2]), x), y), theta[1])+Physics:-diff(diff(diff(Q(x, y, theta[1], theta[2]), x), y), theta[2])-Physics:-diff(Physics:-diff(diff(Q(x, y, theta[1], theta[2]), y), theta[1]), theta[2]) = 0

(1.18)

pde[2] := q[theta[1]] = 0

Physics:-diff(Q(x, y, theta[1], theta[2]), theta[1]) = 0

(1.19)

The solution returned for this system is indeed a general solution

pdsolve([pde[1], pde[2]])

-> Solving ordering for the dependent variables of the PDE system: [_F4(x,y), _F2(x,y), _F3(x,y)]

-> Solving ordering for the independent variables (can be changed using the ivars option): [x, y]
tackling triangularized subsystem with respect to _F4(x,y)
tackling triangularized subsystem with respect to _F2(x,y)
tackling triangularized subsystem with respect to _F3(x,y)
First set of solution methods (general or quasi general solution)
Trying simple case of a single derivative.
HINT = _F6(x)+_F5(y)
Trying HINT = _F6(x)+_F5(y)
HINT is successful
First set of solution methods successful
<- Returning a *general* solution

 

Q(x, y, theta[1], theta[2]) = _F1(x, y)*_lambda1+(_F6(x)+_F5(y))*theta[2]

(1.20)

NULL

This solution involves an anticommutative constant `_&lambda;2`, analogous to the commutative constants _Cn where n is an integer.

 

Download PDE_general_solutions.mw

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


Hi

New developments (after the release of Maple 2016) happened in the project on exact solutions for "Partial Differential Equations & Boundary Conditions". This is work in collaboration with Katherina von Bulow and the improvements are of wide range, representing a noticeable step forward in the capabilities of the Maple system for this kind of problem. As usual, these improvements can be installed in current Maple 2016 by downloading the updated library from the Maplesoft R&D webpage for Differential Equations and Mathematical functions (the update is distributed merged with the updates of the Physics package)

 

The improvements cover:

 

• 

PDE&BC in semi-infinite domains for which a bounded solution is sought

• 

PDE & BC problems in bounded spatial domains via eigenfunction (Fourier) expansions

• 

Implementation of another algebraic method for tackling linear PDE & BC

• 

Improvements in solving PDE & BC solutions by first finding the PDE's general solution.

• 

Improvements in solving PDE & BC problems by using a Fourier transform.

• 

PDE & BC problems that used to require the option HINT = `+` are now solved automatically

 

What follows is a set of examples solved now with these new developments, organized in sections according to the kind of problem. Where relevant, the sections include a subsection on "How it works step by step".

PDE&BC in semi-infinite domains for which a bounded solution is sought can now also be solved via Laplace transforms

 

Maple is now able to solve more PDE&BC problems via Laplace transforms.

 

How it works: Laplace transforms act to change derivatives with respect to one of the independent variables of the domain into multiplication operations in the transformed domain. After applying a Laplace transform to the original problem, we can simplify the problem using the transformed BC, then solve the problem in the transformed domain, and finally apply the inverse Laplace transform to arrive at the final solution. It is important to remember to give pdsolve any necessary restrictions on the variables and constants of the problem, by means of the "assuming" command.

 

A new feature is that we can now tell pdsolve that the dependent variable is bounded, by means of the optional argument HINT = boundedseries.

 

restart

 

Consider the problem of a falling cable lying on a table that is suddenly removed (cf. David J. Logan's Applied Partial Differential Equations p.115).

pde[1] := diff(u(x, t), t, t) = c^2*(diff(u(x, t), x, x))-g
iv[1] := u(x, 0) = 0, u(0, t) = 0, (D[2](u))(x, 0) = 0

 

If we ask pdsolve to solve this problem without the condition of boundedness of the solution, we obtain:

`assuming`([pdsolve([pde[1], iv[1]])], [0 < t, 0 < x, 0 < c])

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

(1.1)

New: If we now ask for a bounded solution, by means of the option HINT = boundedseries, pdsolve simplifies the problem accordingly.

ans[1] := `assuming`([pdsolve([pde[1], iv[1]], HINT = boundedseries)], [0 < t, 0 < x, 0 < c])

u(x, t) = (1/2)*g*(Heaviside(t-x/c)*(c*t-x)^2-c^2*t^2)/c^2

(1.2)

 

And we can check this answer against the original problem, if desired:

`assuming`([pdetest(ans[1], [pde[1], iv[1]])], [0 < t, 0 < x, 0 < c])

[0, 0, 0, 0]

(1.3)

How it works, step by step

 

 Let us see the process this problem undergoes to be solved by pdsolve, step by step.

 

First, the Laplace transform is applied to the PDE:

with(inttrans)

transformed_PDE := laplace((lhs-rhs)(pde[1]), t, s)

s^2*laplace(u(x, t), t, s)-(D[2](u))(x, 0)-s*u(x, 0)-c^2*(diff(diff(laplace(u(x, t), t, s), x), x))+g/s

(1.1.1)

and the result is simplified using the initial conditions:

simplified_transformed_PDE := eval(transformed_PDE, {iv[1]})

s^2*laplace(u(x, t), t, s)-c^2*(diff(diff(laplace(u(x, t), t, s), x), x))+g/s

(1.1.2)

Next, we call the function "laplace(u(x,t),t,s)" by the new name U:

eq_U := subs(laplace(u(x, t), t, s) = U(x, s), simplified_transformed_PDE)

s^2*U(x, s)-c^2*(diff(diff(U(x, s), x), x))+g/s

(1.1.3)

And this equation, which is really an ODE, is solved:

solution_U := dsolve(eq_U, U(x, s))

U(x, s) = exp(-s*x/c)*_F2(s)+exp(s*x/c)*_F1(s)-g/s^3

(1.1.4)

Now, since we want a BOUNDED solution, the term with the positive exponential must be zero, and we are left with:

bounded_solution_U := subs(coeff(rhs(solution_U), exp(s*x/c)) = 0, solution_U)

U(x, s) = exp(-s*x/c)*_F2(s)-g/s^3

(1.1.5)

Now, the initial solution must also be satisfied. Here it is, in the transformed domain:

Laplace_BC := laplace(u(0, t), t, s) = 0

laplace(u(0, t), t, s) = 0

(1.1.6)

Or, in the new variable U,

Laplace_BC_U := U(0, s) = 0

U(0, s) = 0

(1.1.7)

And by applying it to bounded_solution_U, we find the relationship

simplify(subs(x = 0, rhs(bounded_solution_U))) = 0

(_F2(s)*s^3-g)/s^3 = 0

(1.1.8)

isolate((_F2(s)*s^3-g)/s^3 = 0, indets((_F2(s)*s^3-g)/s^3 = 0, unknown)[1])

_F2(s) = g/s^3

(1.1.9)

so that our solution now becomes

bounded_solution_U := subs(_F2(s) = g/s^3, bounded_solution_U)

U(x, s) = exp(-s*x/c)*g/s^3-g/s^3

(1.1.10)

to which we now apply the inverse Laplace transform to obtain the solution to the problem:

`assuming`([u(x, t) = invlaplace(rhs(bounded_solution_U), s, t)], [0 < x, 0 < t, 0 < c])

u(x, t) = (1/2)*g*(-t^2+Heaviside(t-x/c)*(c*t-x)^2/c^2)

(1.1.11)

Four other related examples

 

A few other examples:

pde[2] := diff(u(x, t), t, t) = c^2*(diff(u(x, t), x, x))
iv[2] := u(x, 0) = 0, u(0, t) = g(t), (D[2](u))(x, 0) = 0

ans[2] := `assuming`([pdsolve([pde[2], iv[2]], HINT = boundedseries)], [0 < t, 0 < x, 0 < c])

u(x, t) = Heaviside(t-x/c)*g((c*t-x)/c)

(1.2.1)

`assuming`([pdetest(ans[2], [pde[2], iv[2]])], [0 < t, 0 < x, 0 < c])

[0, 0, 0, 0]

(1.2.2)

pde[3] := diff(u(x, t), t) = k*(diff(u(x, t), x, x)); iv[3] := u(x, 0) = 0, u(0, t) = 1

ans[3] := `assuming`([pdsolve([pde[3], iv[3]], HINT = boundedseries)], [0 < t, 0 < x, 0 < k])

u(x, t) = 1-erf((1/2)*x/(t^(1/2)*k^(1/2)))

(1.2.3)

pdetest(ans[3], [pde[3], iv[3][2]])

[0, 0]

(1.2.4)

pde[4] := diff(u(x, t), t) = k*(diff(u(x, t), x, x)); iv[4] := u(x, 0) = mu, u(0, t) = lambda

ans[4] := `assuming`([pdsolve([pde[4], iv[4]], HINT = boundedseries)], [0 < t, 0 < x, 0 < k])

u(x, t) = (-lambda+mu)*erf((1/2)*x/(t^(1/2)*k^(1/2)))+lambda

(1.2.5)

pdetest(ans[4], [pde[4], iv[4][2]])

[0, 0]

(1.2.6)

 

The following is an example from page 76 in Logan's book:

pde[5] := diff(u(x, t), t) = diff(u(x, t), x, x)
iv[5] := u(x, 0) = 0, u(0, t) = f(t)

ans[5] := `assuming`([pdsolve([pde[5], iv[5]], HINT = boundedseries)], [0 < t, 0 < x])

u(x, t) = (1/2)*x*(int(f(_U1)*exp(-x^2/(4*t-4*_U1))/(t-_U1)^(3/2), _U1 = 0 .. t))/Pi^(1/2)

(1.2.7)

More PDE&BC problems in bounded spatial domains can now be solved via eigenfunction (Fourier) expansions

 

The code for solving PDE&BC problems in bounded spatial domains has been expanded. The method works by separating the variables by product, so that the problem is transformed into an ODE system (with initial and/or boundary conditions) problem, one of which is a Sturm-Liouville problem (a type of eigenvalue problem) which has infinitely many solutions - hence the infinite series representation of the solutions.

restart

 

Here is a simple example for the heat equation:

pde__6 := diff(u(x, t), t) = k*(diff(u(x, t), x, x)); iv__6 := u(0, t) = 0, u(l, t) = 0

ans__6 := `assuming`([pdsolve([pde__6, iv__6])], [0 < l])

u(x, t) = Sum(_C1*sin(_Z1*Pi*x/l)*exp(-k*Pi^2*_Z1^2*t/l^2), _Z1 = 1 .. infinity)

(2.1)

pdetest(ans__6, [pde__6, iv__6])

[0, 0, 0]

(2.2)

 

Now, consider the displacements of a string governed by the wave equation, where c is a constant (cf. Logan p.28).

pde__7 := diff(u(x, t), t, t) = c^2*(diff(u(x, t), x, x))
iv__7 := u(0, t) = 0, u(l, t) = 0

ans__7 := `assuming`([pdsolve([pde__7, iv__7])], [0 < l])

u(x, t) = Sum(sin(_Z2*Pi*x/l)*(sin(c*_Z2*Pi*t/l)*_C1+cos(c*_Z2*Pi*t/l)*_C5), _Z2 = 1 .. infinity)

(2.3)

pdetest(ans__7, [pde__7, iv__7])

[0, 0, 0]

(2.4)

Another wave equation problem (cf. Logan p.130):

pde__8 := diff(u(x, t), t, t)-c^2*(diff(u(x, t), x, x)) = 0; iv__8 := u(0, t) = 0, (D[2](u))(x, 0) = 0, (D[1](u))(l, t) = 0, u(x, 0) = f(x)

ans__8 := `assuming`([pdsolve([pde__8, iv__8], u(x, t))], [0 <= x, x <= l])

u(x, t) = Sum(2*(Int(f(x)*sin((1/2)*Pi*(2*_Z3+1)*x/l), x = 0 .. l))*sin((1/2)*Pi*(2*_Z3+1)*x/l)*cos((1/2)*c*Pi*(2*_Z3+1)*t/l)/l, _Z3 = 1 .. infinity)

(2.5)

pdetest(ans__8, [pde__8, iv__8[1 .. 3]])

[0, 0, 0, 0]

(2.6)

 

Here is a problem with periodic boundary conditions (cf. Logan p.131). The function u(x, t) stands for the concentration of a chemical dissolved in water within a tubular ring of circumference 2*l. The initial concentration is given by f(x), and the variable x is the arc-length parameter that varies from 0 to 2*l.

pde__9 := diff(u(x, t), t) = M*(diff(u(x, t), x, x))
iv__9 := u(0, t) = u(2*l, t), (D[1](u))(0, t) = (D[1](u))(2*l, t), u(x, 0) = f(x)

ans__9 := `assuming`([pdsolve([pde__9, iv__9], u(x, t))], [0 <= x, x <= 2*l])

u(x, t) = (1/2)*_C8+Sum(((Int(f(x)*sin(_Z4*Pi*x/l), x = 0 .. 2*l))*sin(_Z4*Pi*x/l)+(Int(f(x)*cos(_Z4*Pi*x/l), x = 0 .. 2*l))*cos(_Z4*Pi*x/l))*exp(-M*Pi^2*_Z4^2*t/l^2)/l, _Z4 = 1 .. infinity)

(2.7)

pdetest(ans__9, [pde__9, iv__9[1 .. 2]])

[0, 0, 0]

(2.8)

 

The following problem is for heat flow with both boundaries insulated (cf. Logan p.166, 3rd edition)

pde__10 := diff(u(x, t), t) = k*(diff(u(x, t), x, x))
iv__10 := (D[1](u))(0, t) = 0, (D[1](u))(l, t) = 0, u(x, 0) = f(x)

ans__10 := `assuming`([pdsolve([pde__10, iv__10], u(x, t))], [0 <= x, x <= l])

u(x, t) = Sum(2*(Int(f(x)*cos(_Z6*Pi*x/l), x = 0 .. l))*cos(_Z6*Pi*x/l)*exp(-k*Pi^2*_Z6^2*t/l^2)/l, _Z6 = 1 .. infinity)

(2.9)

pdetest(ans__10, [pde__10, iv__10[1 .. 2]])

[0, 0, 0]

(2.10)

 

This is a problem in a bounded domain with the presence of a source. A source term represents an outside influence in the system and leads to an inhomogeneous PDE (cf. Logan p.149):

pde__11 := diff(u(x, t), t, t)-c^2*(diff(u(x, t), x, x)) = p(x, t)
iv__11 := u(0, t) = 0, u(Pi, t) = 0, u(x, 0) = 0, (D[2](u))(x, 0) = 0

ans__11 := pdsolve([pde__11, iv__11], u(x, t))

u(x, t) = Int(Sum(2*(Int(p(x, tau1)*sin(_Z7*x), x = 0 .. Pi))*sin(_Z7*x)*sin(c*_Z7*(t-tau1))/(Pi*_Z7*c), _Z7 = 1 .. infinity), tau1 = 0 .. t)

(2.11)

Current pdetest is unable to verify that this solution cancells the pde__11 mainly because it currently fails in identifying that there is a fourier expansion in it, but its subroutines for testing the boundary conditions work well with this problem

pdetest_BC := `pdetest/BC`

pdetest_BC({ans__11}, [iv__11], [u(x, t)])

[0, 0, 0, 0]

(2.12)

 

 

Consider a heat absorption-radiation problem in the bounded domain 0 <= x and x <= 2, t >= 0:

pde__12 := diff(u(x, t), t) = diff(u(x, t), x, x)
iv__12 := u(x, 0) = f(x), (D[1](u))(0, t)+u(0, t) = 0, (D[1](u))(2, t)+u(2, t) = 0

ans__12 := `assuming`([pdsolve([pde__12, iv__12], u(x, t))], [0 <= x and x <= 2, 0 <= t])

u(x, t) = (1/2)*_C8+Sum(((Int(f(x)*cos((1/2)*_Z8*Pi*x), x = 0 .. 2))*cos((1/2)*_Z8*Pi*x)+(Int(f(x)*sin((1/2)*_Z8*Pi*x), x = 0 .. 2))*sin((1/2)*_Z8*Pi*x))*exp(-(1/4)*Pi^2*_Z8^2*t), _Z8 = 1 .. infinity)

(2.13)

pdetest(ans__12, pde__12)

0

(2.14)

Consider the nonhomogeneous wave equation problem (cf. Logan p.213, 3rd edition):

pde__13 := diff(u(x, t), t, t) = A*x+diff(u(x, t), x, x)
iv__13 := u(0, t) = 0, u(1, t) = 0, u(x, 0) = 0, (D[2](u))(x, 0) = 0

ans__13 := pdsolve([pde__13, iv__13])

u(x, t) = Int(Sum(2*A*(Int(x*sin(Pi*_Z9*x), x = 0 .. 1))*sin(Pi*_Z9*x)*sin(Pi*_Z9*(t-tau1))/(Pi*_Z9), _Z9 = 1 .. infinity), tau1 = 0 .. t)

(2.15)

pdetest_BC({ans__13}, [iv__13], [u(x, t)])

[0, 0, 0, 0]

(2.16)

 

Consider the following Schrödinger equation with zero potential energy (cf. Logan p.30):

pde__14 := I*h*(diff(f(x, t), t)) = -h^2*(diff(f(x, t), x, x))/(2*m)
iv__14 := f(0, t) = 0, f(d, t) = 0

ans__14 := `assuming`([pdsolve([pde__14, iv__14])], [0 < d])

f(x, t) = Sum(_C1*sin(_Z10*Pi*x/d)*exp(-((1/2)*I)*h*Pi^2*_Z10^2*t/(d^2*m)), _Z10 = 1 .. infinity)

(2.17)

pdetest(ans__14, [pde__14, iv__14])

[0, 0, 0]

(2.18)

Another method has been implemented for linear PDE&BC

 

This method is for problems of the form

 

 "(&PartialD;w)/(&PartialD;t)=M[w]"", w(`x__i`,0) = f(`x__i`)" or

 

"((&PartialD;)^2w)/((&PartialD;)^( )t^2)="M[w]", w(`x__i`,0) = f(`x__i`), (&PartialD;w)/(&PartialD;t)() ? ()|() ? (t=0) =g(`x__i`)"

 

where M is an arbitrary linear differential operator of any order which only depends on the spatial variables x__i.

 

Here are some examples:

pde__15 := diff(w(x1, x2, x3, t), t)-(diff(w(x1, x2, x3, t), x2, x1))-(diff(w(x1, x2, x3, t), x3, x1))-(diff(w(x1, x2, x3, t), x3, x3))+diff(w(x1, x2, x3, t), x3, x2) = 0
iv__15 := w(x1, x2, x3, 0) = x1^5*x2*x3NULL

pdsolve([pde__15, iv__15])

w(x1, x2, x3, t) = 20*(((1/20)*x2*x3-(1/20)*t)*x1^2+(1/4)*t*(x2+x3)*x1+t^2)*x1^3

(3.1)

pdetest(%, [pde__15, iv__15])

[0, 0]

(3.2)

 

Here are two examples for which the derivative with respect to t is of the second order, and two initial conditions are given:

pde__16 := diff(w(x1, x2, x3, t), t, t) = diff(w(x1, x2, x3, t), x2, x1)+diff(w(x1, x2, x3, t), x3, x1)+diff(w(x1, x2, x3, t), x3, x3)-(diff(w(x1, x2, x3, t), x3, x2))
iv__16 := w(x1, x2, x3, 0) = x1^3*x2^2+x3, (D[4](w))(x1, x2, x3, 0) = -x2*x3+x1

pdsolve([pde__16, iv__16])

w(x1, x2, x3, t) = x1^3*x2^2+x3-t*x2*x3+t*x1+3*t^2*x2*x1^2+(1/6)*t^3+(1/2)*t^4*x1

(3.3)

pdetest(%, [pde__16, iv__16])

[0, 0, 0]

(3.4)

pde__17 := diff(w(x1, x2, x3, t), t, t) = diff(w(x1, x2, x3, t), x2, x1)+diff(w(x1, x2, x3, t), x3, x1)+diff(w(x1, x2, x3, t), x3, x3)-(diff(w(x1, x2, x3, t), x3, x2))
iv__17 := w(x1, x2, x3, 0) = x1^3*x3^2+sin(x1), (D[4](w))(x1, x2, x3, 0) = cos(x1)-x2*x3

pdsolve([pde__17, iv__17])

w(x1, x2, x3, t) = (1/2)*t^4*x1+t^2*x1^3+3*t^2*x1^2*x3+x1^3*x3^2+(1/6)*t^3-t*x2*x3+cos(x1)*t+sin(x1)

(3.5)

pdetest(%, [pde__17, iv__17])

[0, 0, 0]

(3.6)

More PDE&BC problems are now solved via first finding the PDE's general solution.

 

The following are examples of PDE&BC problems for which pdsolve is successful in first calculating the PDE's general solution, and then fitting the initial or boundary condition to it.

pde__18 := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0
iv__18 := u(0, y) = sin(y)/y

If we ask pdsolve to solve the problem, we get:

ans__18 := pdsolve([pde__18, iv__18])

u(x, y) = (sin(-y+I*x)+_F2(y-I*x)*(y-I*x)+(-y+I*x)*_F2(y+I*x))/(-y+I*x)

(4.1)

and we can check this answer by using pdetest:

pdetest(ans__18, [pde__18, iv__18])

[0, 0]

(4.2)

How it works, step by step:

 

The general solution for just the PDE is:

gensol := pdsolve(pde__18)

u(x, y) = _F1(y-I*x)+_F2(y+I*x)

(4.1.1)

Substituting in the condition iv__18, we get:

u(0, y) = sin(y)/y

(4.1.2)

gensol_with_condition := eval(rhs(gensol), x = 0) = rhs(iv__18)

_F1(y)+_F2(y) = sin(y)/y

(4.1.3)

We then isolate one of the functions above (we can choose either one, in this case), convert it into a function operator, and then apply it to gensol

_F1 = unapply(solve(_F1(y)+_F2(y) = sin(y)/y, _F1(y)), y)

_F1 = (proc (y) options operator, arrow; (-_F2(y)*y+sin(y))/y end proc)

(4.1.4)

eval(gensol, _F1 = (proc (y) options operator, arrow; (-_F2(y)*y+sin(y))/y end proc))

u(x, y) = (-_F2(y-I*x)*(y-I*x)-sin(-y+I*x))/(y-I*x)+_F2(y+I*x)

(4.1.5)

 

 

Three other related examples

 

pde__19 := diff(u(x, y), x, x)+(1/2)*(diff(u(x, y), y, y)) = 0
iv__19 := u(0, y) = sin(y)/y

pdsolve([pde__19, iv__19])

u(x, y) = (2*sin(-y+((1/2)*I)*2^(1/2)*x)+(-I*2^(1/2)*x+2*y)*_F2(y-((1/2)*I)*2^(1/2)*x)+(I*2^(1/2)*x-2*y)*_F2(y+((1/2)*I)*2^(1/2)*x))/(I*2^(1/2)*x-2*y)

(4.2.1)

pdetest(%, [pde__19, iv__19])

[0, 0]

(4.2.2)

pde__20 := diff(u(x, y), x, x)+(1/2)*(diff(u(x, y), y, y)) = 0
iv__20 := u(x, 0) = sin(x)/x

pdsolve([pde__20, iv__20])

u(x, y) = (sinh((1/2)*(I*2^(1/2)*x-2*y)*2^(1/2))*2^(1/2)-(I*2^(1/2)*x-2*y)*(_F2(-y+((1/2)*I)*2^(1/2)*x)-_F2(y+((1/2)*I)*2^(1/2)*x)))/(I*2^(1/2)*x-2*y)

(4.2.3)

pdetest(%, [pde__20, iv__20])

[0, 0]

(4.2.4)

pde__21 := diff(u(r, t), r, r)+(diff(u(r, t), r))/r+(diff(u(r, t), t, t))/r^2 = 0
iv__21 := u(3, t) = sin(6*t)

ans__21 := pdsolve([pde__21, iv__21])

u(r, t) = -_F2(-(2*I)*ln(3)+I*ln(r)+t)+sin(-(6*I)*ln(3)+(6*I)*ln(r)+6*t)+_F2(-I*ln(r)+t)

(4.2.5)

pdetest(ans__21, [pde__21, iv__21])

[0, 0]

(4.2.6)

More PDE&BC problems are now solved by using a Fourier transform.

 

restart

Consider the following problem with an initial condition:

pde__22 := diff(u(x, t), t) = diff(u(x, t), x, x)+m
iv__22 := u(x, 0) = sin(x)

 

pdsolve can solve this problem directly:

ans__22 := pdsolve([pde__22, iv__22])

u(x, t) = sin(x)*exp(-t)+m*t

(5.1)

And we can check this answer against the original problem, if desired:

pdetest(ans__22, [pde__22, iv__22])

[0, 0]

(5.2)

How it works, step by step

 

Similarly to the Laplace transform method, we start the solution process by first applying the Fourier transform to the PDE:

with(inttrans)

transformed_PDE := fourier((lhs-rhs)(pde__22) = 0, x, s)

-2*m*Pi*Dirac(s)+s^2*fourier(u(x, t), x, s)+diff(fourier(u(x, t), x, s), t) = 0

(5.1.1)

Next, we call the function "fourier(u(x,t),x,s1)" by the new name U:

transformed_PDE_U := subs(fourier(u(x, t), x, s) = U(t, s), transformed_PDE)

-2*m*Pi*Dirac(s)+s^2*U(t, s)+diff(U(t, s), t) = 0

(5.1.2)

And this equation, which is really an ODE, is solved:

solution_U := dsolve(transformed_PDE_U, U(t, s))

U(t, s) = (2*m*Pi*Dirac(s)*t+_F1(s))*exp(-s^2*t)

(5.1.3)

Now, we apply the Fourier transform to the initial condition iv__22:

u(x, 0) = sin(x)

(5.1.4)

transformed_IC := fourier(iv__22, x, s)

fourier(u(x, 0), x, s) = I*Pi*(Dirac(s+1)-Dirac(s-1))

(5.1.5)

Or, in the new variable U,

trasnformed_IC_U := U(0, s) = rhs(transformed_IC)

U(0, s) = I*Pi*(Dirac(s+1)-Dirac(s-1))

(5.1.6)

Now, we evaluate solution_U at t = 0:

solution_U_at_IC := eval(solution_U, t = 0)

U(0, s) = _F1(s)

(5.1.7)

and substitute the transformed initial condition into it:

eval(solution_U_at_IC, {trasnformed_IC_U})

I*Pi*(Dirac(s+1)-Dirac(s-1)) = _F1(s)

(5.1.8)

Putting this into our solution_U, we get

eval(solution_U, {(rhs = lhs)(I*Pi*(Dirac(s+1)-Dirac(s-1)) = _F1(s))})

U(t, s) = (2*m*Pi*Dirac(s)*t+I*Pi*(Dirac(s+1)-Dirac(s-1)))*exp(-s^2*t)

(5.1.9)

Finally, we apply the inverse Fourier transformation to this,

solution := u(x, t) = invfourier(rhs(U(t, s) = (2*m*Pi*Dirac(s)*t+I*Pi*(Dirac(s+1)-Dirac(s-1)))*exp(-s^2*t)), s, x)

u(x, t) = sin(x)*exp(-t)+m*t

(5.1.10)

PDE&BC problems that used to require the option HINT = `+` to be solved are now solved automatically

 

The following two PDE&BC problems used to require the option HINT = `+` in order to be solved. This is now done automatically within pdsolve.

pde__23 := diff(u(r, t), r, r)+(diff(u(r, t), r))/r+(diff(u(r, t), t, t))/r^2 = 0
iv__23 := u(1, t) = 0, u(2, t) = 5

ans__23 := pdsolve([pde__23, iv__23])

u(r, t) = 5*ln(r)/ln(2)

(6.1)

pdetest(ans__23, [pde__23, iv__23])

[0, 0, 0]

(6.2)

pde__24 := diff(u(x, y), y, y)+diff(u(x, y), x, x) = 6*x-6*y

iv__24 := u(x, 0) = x^3+11*x+1, u(x, 2) = x^3+11*x-7, u(0, y) = -y^3+1, u(4, y) = -y^3+109

ans__24 := pdsolve([pde__24, iv__24])

u(x, y) = x^3-y^3+11*x+1

(6.3)

pdetest(ans__24, [pde__24, iv__24])

[0, 0, 0, 0, 0]

(6.4)

``



Download PDE_and_BC_update.mw

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

Greetings to all.

With the following matter I am betting on there being a simple mistake on my part due to fatigue owing to a challenging session of intense computing. The following link at Math.Stackexchange.Com points to a computation involving complex residues. Consult the link for additional details.

I usually verify my computations with Maple, I did the same this time. Thereby I happened on a curious phenomenon which I have documented below. Please study the session data provided, I believe it speaks for itself.

user@host:~$ math
Mathematica 10.0 for Linux x86 (64-bit)
Copyright 1988-2014 Wolfram Research, Inc.

In[1]:= Residue[z^2/(z^4 + 2*z^2 + 2)^2, {z, 2^(1/4)*Exp[3*Pi*I/8]}]

            1/8      1/4
        (-1)    ((-1)    + Sqrt[2])
Out[1]= ----------------------------
            1/4      1/4           3
        16 2    ((-1)    - Sqrt[2])

In[2]:= N[%]

Out[2]= 0.117223 - 0.0083308 I

In[3]:=
user@host:~$ maple
    |\^/|     Maple 18 (X86 64 LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2014
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> residue(z^2/(z^4 + 2*z^2 + 2)^2, z=2^(1/4)*exp(3*Pi*I/8));
                                       0

> quit
memory used=0.9MB, alloc=8.3MB, time=0.07
user@host:~$

I am looking forward to learning what the correct syntax is to get the residue in this case and I hope I can assist other users who might have run into the same problem. I will cancel the question should it turn out to be trivial and of little potential use to the community.

Best regards,

Marko Riedel

Post Scriptum. Being a programmer myself I would be curious to learn more about the algorithms that are deployed here and how and why they did not succeed.

Formulating and solving the equivalence problem for Schwarzschild metric in a simple case

 

In connection with the digitizing in Maple 2016 of the database of solutions to Einstein's equations of the book Exact Solutions to Einstein Field Equations. I was recently asked about a statement found in the "What is new in Physics in Maple 2016" page:

  

In the Maple PDEtools package, you have the mathematical tools - including a complete symmetry approach - to work with the underlying [Einstein’s] partial differential equations. [By combining that functionality with the one in the Physics and Physics:-Tetrads package] you can also formulate and, depending on the metrics also resolve, the equivalence problem; that is: to answer whether or not, given two metrics, they can be obtained from each other by a transformation of coordinates, as well as compute the transformation.

This question posed is a reasonable one: "could you please provide one example?" This post provides that example.

 

First of all the existing science behind: in my opinion, the main reference regarding the equivalence problem is at the paper "A Review of the Geometrical Equivalence of Metrics in General Relativity", General Relativity and Gravitation, Vol. 12, No. 9, 1980, by A. Karlhede (University of Stockholm). This approach got refined later by others and, generally speaking, it is currently know as the Cartan-Karlhede method, summarized in chapter 9.2 of the book Exact Solutions to Einstein Field Equations. whose solutions were all digitized within the Physics and DifferentialGeometry packages for Maple 2016. This method of Chapter 9.2 (see also Tetrads and Weyl scalars in canonical form, Mapleprimes post), however, is not the only approach to the problem, and sometimes simpler methods can handle the problem faster, or just in simpler forms.

 

The example worked out below is actually the example from Karlhede's paper just mentioned, on pages 704 - 706: "Show that the Schwarzschild metric and its form written in terms of isotropic spherical coordinates are equivalent, and derive the transformation that relates them". Because this problem happens to be simple for nowadays computer algebra, below I also tackle it modified, slightly more difficult variants of it. The approach shown works for more complicated cases as well.

 

Below we tackle Karlhede's paper-problem using: one PDEtools command, the Physics:-TransformCoordinates, the Physics:-Weyl command to compute the Weyl scalars and the Physics:-Tetrads:-PetrovType to see the Petrov type of the metrics involved. The transformation resolving the equivalence is explicitly derived.

 

Start loading the Physics and Tetrads package. To reproduce the computations below, as usual, update your Physics library with the one available for download at the Maplesoft R&D Physics webpage

with(Physics); with(Tetrads); Setup(auto = true, tetradmetric = null, signature = `+---`)

`Setting lowercaselatin letters to represent tetrad indices `

 

0, "%1 is not a command in the %2 package", Tetrads, Physics

 

0, "%1 is not a command in the %2 package", Tetrads, Physics

 

`* Partial match of  'auto' against keyword 'automaticsimplification'`

 

[automaticsimplification = true, signature = `+ - - -`, tetradmetric = {(1, 2) = 1, (3, 4) = -1}]

(1)

To formulate the problem, set first some symbols to represent the changed metric, changed mass and changed coordinates - no mathematics at this point

gt, mt, tt, rt, thetat, phit := `&gfr;`, `&mfr;`, `&tfr;`, `&rfr;`, `&vartheta;`, `&varphi;`

`&gfr;`, `&mfr;`, `&tfr;`, `&rfr;`, vartheta, varphi

(2)

Set now a new coordinates system, call it Y, involving the new coordinates (in the paper they are represented with a tilde on top of the letters)

Coordinates(Y = [tt, rt, thetat, phit])

`Default differentiation variables for d_, D_ and dAlembertian are: `*{Y = (`&tfr;`, `&rfr;`, `&vartheta;`, `&varphi;`)}

 

`Systems of spacetime Coordinates are: `*{Y = (`&tfr;`, `&rfr;`, `&vartheta;`, `&varphi;`)}

 

{Y}

(3)

According to eq.(7.6) of the paper, the line element of Schwarzschild solution in isotropic spherical coordinates is given by

`#msup(mi("ds"),mn("2"))` := ((1-mt/(2*rt))/(1+mt/(2*rt)))^2*d_(tt)^2-(1+mt/(2*rt))^4*(d_(rt)^2+rt^2*d_(thetat)^2+rt^2*sin(thetat)^2*d_(phit)^2)

(-2*`&rfr;`+`&mfr;`)^2*Physics:-d_(`&tfr;`)^2/(2*`&rfr;`+`&mfr;`)^2-(1/16)*(2*`&rfr;`+`&mfr;`)^4*(Physics:-d_(`&rfr;`)^2+`&rfr;`^2*Physics:-d_(vartheta)^2+`&rfr;`^2*sin(vartheta)^2*Physics:-d_(varphi)^2)/`&rfr;`^4

(4)

Set this to be the metric

Setup(metric = `#msup(mi("ds"),mn("2"))`)

Check it out

g_[]

Physics:-g_[mu, nu] = Matrix(%id = 18446744078306516254)

(5)

In connection with the transformation used further below, compute now the Petrov type and the Weyl scalars for this metric, just to have an idea of what is behind this metric.

PetrovType()

"D"

(6)

Weyl[scalars]

psi__0 = 0, psi__1 = 0, psi__2 = -64*`&rfr;`^3*`&mfr;`/(2*`&rfr;`+`&mfr;`)^6, psi__3 = 0, psi__4 = 0

(7)

We see that the Weyl scalars are already in canonical form (see post in Mapleprimes about canonical forms): only `&Psi;__2` <> 0 and the important thing: it depends on only one coordinate, `&rfr;` .

 

Now: we want to see if this metric (5) is equivalent to Schwarzschild metric in standard spherical coordinates

g_[sc]

`Systems of spacetime Coordinates are: `*{X = (t, r, theta, phi), Y = (`&tfr;`, `&rfr;`, `&vartheta;`, `&varphi;`)}

 

`Default differentiation variables for d_, D_ and dAlembertian are: `*{X = (t, r, theta, phi)}

 

`The Schwarzschild metric in coordinates `[t, r, theta, phi]

 

`Parameters: `[m]

 

Physics:-g_[mu, nu] = Matrix(%id = 18446744078795590102)

(8)

The equivalence we want to resolve is regarding an arbitrary relationship `&mfr;`(m)between the masses used in (5) and (8) and a generic change of variables from X to Y

TR := {phi = Phi(Y), r = R(Y), t = Tau(Y), theta = Theta(Y)}

{phi = Phi(Y), r = R(Y), t = Tau(Y), theta = Theta(Y)}

(9)

Using a differential equation mindset, the formulation of the equivalence between (8) and (5) under the transformation (9) is actually simple: change variables in (8), using (9) and the Physics:-TransformCoordinates command (this is the command that changes variables in tensorial expressions), then equate the result to (5), then try to solve the problem for the unknowns `&mfr;`(m), Phi(Y), R(Y), Theta(Y) and Tau(Y).

 

We note at this point, however, that the Weyl scalars for Schwarzschild metric in this standard form (8) are also in canonical form of Petrov type D and also depend on only one variable, r 

PetrovType()

"D"

(10)

Weyl[scalars]

psi__0 = 0, psi__1 = 0, psi__2 = -m/r^3, psi__3 = 0, psi__4 = 0

(11)

The fact that the Weyl scalars in both cases ((7) and (11)) are in canonical form (only `&Psi;__2` <> 0 ) and in both cases this scalar depends on only one coordinate is already an indicator that the transformation involved changes only one variable in terms of the other one. So one could just search for a transformation of the form r = R(`&rfr;`) and resolve the problem instantly. Still, to make the problem slightly more general, consider instead a generic transformation for r in terms of all of Y = (`&tfr;`, `&rfr;`, `&vartheta;`, `&varphi;`)

tr := r = R(Y)

r = R(Y)

(12)

PDEtools:-declare(r = R(Y))

R(`&tfr;`, `&rfr;`, vartheta, varphi)*`will now be displayed as`*R

(13)

Transform the  coordinates in the metric (because of having used PDEtools:-declare, derivatives of the unknowns R are displayed indexed, for compact notation)

TransformCoordinates(tr, g_[mu, nu])

Matrix(%id = 18446744078873927542)

(14)

Proceed equating (14) to (5) to obtain a set of equations that entirely formulates the problem

"convert(rhs(?)=? ,setofequations)"

{0 = (diff(R(Y), `&rfr;`))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), vartheta))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), (-2*`&rfr;`+`&mfr;`)^2/(2*`&rfr;`+`&mfr;`)^2 = ((diff(R(Y), `&tfr;`))^2*R(Y)^2-4*(-(1/2)*R(Y)+m)^2)/(R(Y)*(-R(Y)+2*m)), -(1/16)*(2*`&rfr;`+`&mfr;`)^4/`&rfr;`^4 = (diff(R(Y), `&rfr;`))^2*R(Y)/(-R(Y)+2*m), -(1/16)*(2*`&rfr;`+`&mfr;`)^4/`&rfr;`^2 = -(diff(R(Y), vartheta))^2*R(Y)/(R(Y)-2*m)-R(Y)^2, -(1/16)*(2*`&rfr;`+`&mfr;`)^4*sin(vartheta)^2/`&rfr;`^2 = 2*((1/2)*(diff(R(Y), varphi))^2+(cos(vartheta)-1)*R(Y)*(cos(vartheta)+1)*(-(1/2)*R(Y)+m))*R(Y)/(-R(Y)+2*m)}

(15)

This problem, shown in Karlhede's paper as the example of the approach he summarized, is solvable using the differential equation commands of PDEtools (in this case casesplit) in one go and no time, obtaining the same solution shown in the paper with equation number (7.10), the problem actually admits two solutions

PDEtools:-casesplit({0 = (diff(R(Y), `&rfr;`))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), vartheta))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), (-2*`&rfr;`+`&mfr;`)^2/(2*`&rfr;`+`&mfr;`)^2 = ((diff(R(Y), `&tfr;`))^2*R(Y)^2-4*(-(1/2)*R(Y)+m)^2)/(R(Y)*(-R(Y)+2*m)), -(1/16)*(2*`&rfr;`+`&mfr;`)^4/`&rfr;`^4 = (diff(R(Y), `&rfr;`))^2*R(Y)/(-R(Y)+2*m), -(1/16)*(2*`&rfr;`+`&mfr;`)^4/`&rfr;`^2 = -(diff(R(Y), vartheta))^2*R(Y)/(R(Y)-2*m)-R(Y)^2, -(1/16)*(2*`&rfr;`+`&mfr;`)^4*sin(vartheta)^2/`&rfr;`^2 = 2*((1/2)*(diff(R(Y), varphi))^2+(cos(vartheta)-1)*R(Y)*(cos(vartheta)+1)*(-(1/2)*R(Y)+m))*R(Y)/(-R(Y)+2*m)}, [R, mt])

`casesplit/ans`([R(Y) = -(1/4)*(m-2*`&rfr;`)^2/`&rfr;`, `&mfr;` = -m], []), `casesplit/ans`([R(Y) = (1/4)*(2*`&rfr;`+m)^2/`&rfr;`, `&mfr;` = m], [])

(16)

By all means this does not mean this differential equation approach is better than the general approach mentioned in the paper (also in section 9.2 of the Exact Solutions book). This presentation above only makes the point of the paragraph mentioned at the beginning of this worksheet "... [in Maple 2016] you can also formulate and, depending on the the metrics also resolve, the equivalence problem; that is: to answer whether or not, given two metrics, they can be obtained from each other by a transformation of coordinates, as well as compute the transformation." 

 

In any case this problem above is rather easy for the computer. Consider a slightly more difficult problem, where `&mfr;` <> m. For example:

"subs(mt = 1/(mt^(2)),?)"

Physics:-g_[mu, nu] = Matrix(%id = 18446744078854733566)

(17)

Tackle now the same problem

"convert(rhs(?)=? ,setofequations)"

{0 = (diff(R(Y), `&rfr;`))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), vartheta))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), (-2*`&rfr;`+1/`&mfr;`^2)^2/(2*`&rfr;`+1/`&mfr;`^2)^2 = ((diff(R(Y), `&tfr;`))^2*R(Y)^2-4*(-(1/2)*R(Y)+m)^2)/(R(Y)*(-R(Y)+2*m)), -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4/`&rfr;`^4 = (diff(R(Y), `&rfr;`))^2*R(Y)/(-R(Y)+2*m), -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4/`&rfr;`^2 = -(diff(R(Y), vartheta))^2*R(Y)/(R(Y)-2*m)-R(Y)^2, -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4*sin(vartheta)^2/`&rfr;`^2 = 2*((1/2)*(diff(R(Y), varphi))^2+(cos(vartheta)-1)*R(Y)*(cos(vartheta)+1)*(-(1/2)*R(Y)+m))*R(Y)/(-R(Y)+2*m)}

(18)

The solutions to the equivalence between (17) and (5) are then given by

PDEtools:-casesplit({0 = (diff(R(Y), `&rfr;`))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), varphi))*(diff(R(Y), vartheta))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&rfr;`))*R(Y)/(-R(Y)+2*m), 0 = (diff(R(Y), vartheta))*(diff(R(Y), `&tfr;`))*R(Y)/(-R(Y)+2*m), (-2*`&rfr;`+1/`&mfr;`^2)^2/(2*`&rfr;`+1/`&mfr;`^2)^2 = ((diff(R(Y), `&tfr;`))^2*R(Y)^2-4*(-(1/2)*R(Y)+m)^2)/(R(Y)*(-R(Y)+2*m)), -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4/`&rfr;`^4 = (diff(R(Y), `&rfr;`))^2*R(Y)/(-R(Y)+2*m), -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4/`&rfr;`^2 = -(diff(R(Y), vartheta))^2*R(Y)/(R(Y)-2*m)-R(Y)^2, -(1/16)*(2*`&rfr;`+1/`&mfr;`^2)^4*sin(vartheta)^2/`&rfr;`^2 = 2*((1/2)*(diff(R(Y), varphi))^2+(cos(vartheta)-1)*R(Y)*(cos(vartheta)+1)*(-(1/2)*R(Y)+m))*R(Y)/(-R(Y)+2*m)}, [R, mt])

`casesplit/ans`([R(Y) = -(1/4)*(m-2*`&rfr;`)^2/`&rfr;`, `&mfr;`^2 = -1/m], [`&mfr;` <> 0]), `casesplit/ans`([R(Y) = (1/4)*(2*`&rfr;`+m)^2/`&rfr;`, `&mfr;`^2 = 1/m], [`&mfr;` <> 0])

(19)

Moreover, despite that the Weyl scalars suggest that a transformation of only one variable is sufficient to solve the problem, one could also consider a more general transformation, of more variables. Provided we exclude theta (because there is cos(theta) around and that would take us to solve differential equations for Theta(theta), that involve things like cos(Theta(theta))), and also to speed up matters let's remove the change in phi, consider an arbitrary change in r and t

TR := select(has, {phi = Phi(Y), r = R(Y), t = Tau(Y), theta = Theta(Y)}, {r, t})

{r = R(Y), t = Tau(Y)}

(20)

PDEtools:-declare({r = R(Y), t = Tau(Y)})

R(`&tfr;`, `&rfr;`, vartheta, varphi)*`will now be displayed as`*R