Solving a system of linear PDEs using Mathematica

Question

I have the following system of PDEs:

Where: phi(d,e,f,g,h,j,k)

and c1-c8 are real valued constants

I coded it into mathematica like this:

eqns = {C1*k*D[phi[d, e, f, g, h, j, k], g] - 
    C2*f*D[phi[d, e, f, g, h, j, k], k] + 
    C1*k*D[phi[d, e, f, g, h, j, k], f] - 
    C3*k*D[phi[d, e, f, g, h, j, k], d] == 0,
      C1*k*D[phi[d, e, f, g, h, j, k], h] - 
    C2*g*D[phi[d, e, f, g, h, j, k], k] + 
    C4*k*D[phi[d, e, f, g, h, j, k], j] - 
    C3*k*D[phi[d, e, f, g, h, j, k], e] == 0,
      C5*k*D[phi[d, e, f, g, h, j, k], d] - 
    C2*h*D[phi[d, e, f, g, h, j, k], k] - 
    C6*k*D[phi[d, e, f, g, h, j, k], g] - 
    C6*k*D[phi[d, e, f, g, h, j, k], f] == 0,
      C5*k*D[phi[d, e, f, g, h, j, k], e] - 
    C2*j*D[phi[d, e, f, g, h, j, k], k] - 
    C7*k*D[phi[d, e, f, g, h, j, k], j] + 
    C8*k*D[phi[d, e, f, g, h, j, k], j] - 
    C6*k*D[phi[d, e, f, g, h, j, k], h] == 0}

And called Dsolve:

Dsolve[eqns, phi[d, e, f, g, h, j, k], phi[d, e, f, g, h, j, k]]

Unfortunately Dsolve is not able to solve it.

I formulated these PDEs from an engineering problem and to the best of my knowledge it should have atleast 4 solutions to it.

I want to know if this is beyond mathematica's ability to solve or if I did some preprocesssing of some kind if it would help?

I am just starting to learn about PDEs so I'm not quite sure what I could do from my end to help Mathematica solve it.

Answer 1

DSolve has difficulty solving systems of first-order linear PDEs, for instance, 163482 or 174971. Here, too it can obtain a solution, but not without assistance. Begin by solving the first of the four equations.

s1 = DSolveValue[eqns[[1]], phi[d, e, f, g, h, j, k], {d, e, f, g, h, j, k}] /. 
    C[1][z1__][z2__] -> xi1[z1, z2]
(* xi1[e, h, j, (C1 d + C3 f)/C3, (C1 d + C3 g)/C3, 
   (-C1 C2 d^2 - 2 C2 C3 d f + C3^2 k^2)/(2 C3^2)] *)

Then insert this result into the fourth equation.

Simplify[#/k & /@ (eqns[[4]] /. phi -> Function[{d, e, f, g, h, j, k}, Evaluate@s1])]
(* C2*j*Derivative[0, 0, 0, 0, 0, 1][xi1][e, h, j, (C1*d)/C3 + f, (C1*d)/C3 + g, 
   -(C1*C2*d^2 + 2*C2*C3*d*f - C3^2*k^2)/(2*C3^2)] + ... *)

(The whole expression is a bit long to reproduce here.) Now, DSolve does not like independent variable like (C1*d)/C3 + f, so they must be replaced by dummy variables.

var = Array[Unique["v"] &, 6];
rls = Simplify[Thread[List @@ s1 -> var]]
(* {e -> v11, h -> v12, j -> v13, (C1 d)/C3 + f -> v14, (C1 d)/C3 + g -> v15, 
   -((C1 C2 d^2 + 2 C2 C3 d f - C3^2 k^2)/(2 C3^2)) -> v16} *)

s4 = Simplify[DSolveValue[%%% /. rls, xi1 @@ var, var] /. 
    C[1][z1__][z2__] -> xi4[z1, z2] /. (Reverse[#] & /@ rls)]

(* xi4[(C1 d)/C3 + f, (C1 d)/C3 + g, (C6 e)/C5 + h, (C7 e - C8 e + C5 j)/C5, 
   1/2 (-((C1 C2 d^2)/C3^2) - (2 C2 d f)/C3 + (C2 e (C7 e - C8 e + 2 C5 j))/C5^2
   + k^2)] *)

Substituting this result into the third equation then gives

Simplify[#/k & /@ (eqns[[3]] /. phi -> Function[{d, e, f, g, h, j, k}, Evaluate@s4])]
(* (C2*(C1*C5*d + C3*(-(C6*d) + C5*f + C3*h))*Derivative[0, 0, 0, 0, 1][xi4]
   [(C1*d)/C3 + f, (C1*d)/C3 + g, (C6*e)/C5 + h, (C7*e - C8*e + C5*j)/C5, 
   (-((C1*C2*d^2)/C3^2) - (2*C2*d*f)/C3 + (C2*e*(C7*e - C8*e + 2*C5*j))/C5^2
   + k^2)/2])/C3 + ... *)

To proceed, the coefficient (C2*(C1*C5*d + C3*(-(C6*d) + C5*f + C3*h)) must be expressed in terms of the six compound arguments of xi4 or eliminated some other way. In fact, it is easy to see that the coefficient cannot be expressed in terms of the six arguments, so instead the derivative must be set equal to zero. This is accomplished mechanically by

% /. Derivative[n__, 1][xi4][z__, _] -> 0 /. 
    Derivative[n__, 0][xi4][z__, _] -> Derivative[n][xi4][z]
(* (-(C1*C5) + C3*C6)*(Derivative[0, 1, 0, 0][xi4][(C1*d)/C3 + f, (C1*d)/C3 + g, 
   (C6*e)/C5 + h, (C7*e - C8*e + C5*j)/C5] + Derivative[1, 0, 0, 0][xi4]
   [(C1*d)/C3 + f, (C1*d)/C3 + g, (C6*e)/C5 + h, (C7*e - C8*e + C5*j)/C5]) == 0 *)

Now, solve as with the earlier intermediate results.

var = Array[Unique["v"] &, 4];
rls = Simplify[Thread[Most[List @@ s4] -> var]];
s3 = Simplify[DSolveValue[%%% /. rls, xi4 @@ var, var] /. 
    C[1][z1__][z2__] -> xi3[z1, z2] /. (Reverse[#] & /@ rls)]
(* xi3[(C6 e)/C5 + h, (C7 e - C8 e + C5 j)/C5, -f + g] *)

Finally, insert this into the remaining equation and solve.

Simplify[eqns[[2]] /. phi -> Function[{d, e, f, g, h, j, k}, Evaluate@s3]];
var = Array[Unique["v"] &, 3];
rls = Simplify[Thread[List @@ s3 -> var]]
s2 = Simplify[DSolveValue[%%% /. rls, xi3 @@ var, var] /. 
    C[1][z1__][z2__] -> xi2[z1, z2] /. (Reverse[#] & /@ rls)]
(* xi2[-f + g, -((C4 C6 e - C1 C7 e + C1 C8 e + C4 C5 h - C3 C7 h + C3 C8 h - 
   C1 C5 j + C3 C6 j)/(C1 C5 - C3 C6))] *)

which is the desired solution. phi is an arbitrary function of the two remaining compound variables. To verify this result, substitute it into eqns.

Simplify[eqns /. phi -> Function[{d, e, f, g, h, j, k}, Evaluate@s2]]
(* {True, True, True, True} *)

Answer 2

The answer to question 32843 provides an interesting alternative approach for solving the question here. Beginning with the left sides of the equations in the question,

eqn = First /@ eqns;

recast the expressions into an equivalent form that can be integrated readily.

eqn2 = Union[eqn, Flatten[D[eqn, #] & /@ {d, e, f, g, h, j, k}]] // Simplify;
vars = Select[Variables[eqn2], ! FreeQ[#, phi] &];
elim = Cases[vars, Derivative[z__][_][__] /; Total[{z}] > 1];
keep = Complement[vars, elim];
gb = GroebnerBasis[Flatten@eqn2, keep, elim, CoefficientDomain -> 
    RationalFunctions, MonomialOrder -> EliminationOrder]

(* {Derivative[1, 0, 0, 0, 0, 0, 0][phi][d, e, f, g, h, j, k], 
    Derivative[0, 0, 0, 1, 0, 0, 0][phi][d, e, f, g, h, j, k] + 
      Derivative[0, 0, 1, 0, 0, 0, 0][phi][d, e, f, g, h, j, k], 
    (-(C4*C6) + C1*C7 - C1*C8)*Derivative[0, 0, 0, 0, 1, 0, 0][phi][d, e, f, g, h, j, k] + 
      (C4*C5 - C3*C7 + C3*C8)*Derivative[0, 1, 0, 0, 0, 0, 0][phi][d, e, f, g, h, j, k], 
    (C4*C6 - C1*C7 + C1*C8)*Derivative[0, 0, 0, 0, 0, 1, 0][phi][d, e, f, g, h, j, k] + 
      (C1*C5 - C3*C6)*Derivative[0, 1, 0, 0, 0, 0, 0][phi][d, e, f, g, h, j, k], 
    Derivative[0, 0, 0, 0, 0, 0, 1][phi][d, e, f, g, h, j, k]} *)

By inspection, phy must be an arbitrary function of g - f (from gb[[2]]) and of some linear combination of e, g, h (from gb[[3;;4]]), and independent of d (from gb[[1]]) and k (from gb[[5]]). The linear combination can be determined by

Coefficient[gb /. phi -> 
Function[{d, e, f, g, h, j, k}, x4[ae e + ah h + aj j, g - f]], 
    Derivative[1, 0][x4][ae*e + ah*h + aj*j, -f + g]
Flatten@Solve[Thread[% == 0], {ae, ah, aj}]

(* {ah -> -((ae (-C4 C5 + C3 C7 - C3 C8))/(C4 C6 - C1 C7 + C1 C8)), 
    aj -> -((ae (-C1 C5 + C3 C6))/(-C4 C6 + C1 C7 - C1 C8))} *)

This is equivalent to my earlier answer, obtained with greater effort.

If desired, this answer also can be obtained by systematically solving the five PDEs embodied by gb. First solve gb[[1;;5;;4]].

gb1 = DeleteCases[gb /. (Flatten@DSolve[Thread[Drop[gb, {2, 4}] == 0], 
    phi, {d, e, f, g, h, j, k}] /. C[1] -> x1), 0];

Next, solve gb1[[1]], equivalent to gb[[2]].

gb2 = DeleteCases[gb1 /. (Flatten@DSolve[First@gb1 == 0, x1, {e, f, g, h, j}] 
    /. C[1][z1__][z2_] -> x2[z1, z2]), 0]

(* {(-(C4*C6) + C1*C7 - C1*C8)*Derivative[0, 1, 0, 0][x2][e, h, j, -f + g] + 
      (C4*C5 - C3*C7 + C3*C8)*Derivative[1, 0, 0, 0][x2][e, h, j, -f + g], 
    (C4*C6 - C1*C7 + C1*C8)*Derivative[0, 0, 1, 0][x2][e, h, j, -f + g] + 
      (C1*C5 - C3*C6)*Derivative[1, 0, 0, 0][x2][e, h, j, -f + g]} *)

and solve these remaining two coupled PDEs in three variables by the method employed in my earlier answer or the method employed earlier in this answer. Incidentally, ae*e + ah*h + aj*j and -f + g can be viewed as characteristics of the original PDEs.