# Solving a system of partial differential equations consist of 6 equations on 9 variables by using Mathematica

I have a system of partial differential equations consist of 6 equations on 9 variables.

p = D[f[x1, x2, x3, y1, y2, y3, z1, z2, z3], x1];

q = D[f[x1, x2, x3, y1, y2, y3, z1, z2, z3], x2];

r = D[f[x1, x2, x3, y1, y2, y3, z1, z2, z3], x3];

o = D[f[x1, x2, x3, y1, y2, y3, z1, z2, z3], y1];

x = D[f[x1, x2, x3, y1, y2, y3, z1, z2, z3], y2];

a = D[f[x1, x2, x3, y1, y2, y3, z1, z2, z3], y3];

b = D[f[x1, x2, x3, y1, y2, y3, z1, z2, z3], z1];

c = D[f[x1, x2, x3, y1, y2, y3, z1, z2, z3], z2];

d = D[f[x1, x2, x3, y1, y2, y3, z1, z2, z3], z3];

equ1 = 2 x1 p + 2 x2 q + 2 x3 r - y1 o - y2 x - y3 a;

equ2 = -x1 p - x2 q - x3 r - z1 b - z2 c - z3 d + 2 y1 o + 2 y2 x + 2 y3 a;

equ3 = 2 z1 b + 2 z2 c + 2 z3 d - y1 o - y2 x - y3 a;

equ4 = (x1 (x1 + 2 x2 + 2 x3)) p + (x2 (x2 +
2 x3)) q + (x3^2) r - (y1 (x2 + x3)) o - y2 x3 x;

equ5 = (y1 (y1 + 2 y2 + 2 y3)) o + (y2 (y2 + 2 y3)) x + (y3^2) a - (x1 (y1 + y2 + y3)) p
- x2 (y2 + y3) q - x3 y3 r  - z1 (y2 + y3) b - z2 y3 c;

equ6 = (z1 (z1 + 2 z2 + 2 z3)) b + (z2 (z2 + 2 z3)) c + (z3^2) d - (y1 (z1 + z2 + z3)) o
- y2 ( z2 + z3) x - y3 z3 a;


The first three equations are degree equations and show that the solution has degree zero. so somehow by changing the variables to

v1=x2 /x1, v2=x3 /x1, w1=y2/y1, w2=y3/y1, k1=z2/z1, k2=z3/z1,


we will have 3 equations with 6 variables.

Can someone give me a more clever idea than using Dsolve for such kind of systems? because Dsolve returns un-evaluated.

DSolve[{equ1 == 0, equ2 == 0, equ3 == 0, equ4 == 0, equ5 == 0, equ6 == 0},
f[x1, x2, x3, y1, y2, y3, z1, z2, z3], {x1, x2, x3, y1, y2, y3, z1, z2, z3}]


Thanks!

• Why do you think you can find nine variables from six equations, or six variables from three equations? What about two variables from one equation (e.g., $x + y = 7$)? Dec 16, 2015 at 16:55
• In principle being underdetermined there could be many solutions. There is really no chance you will find an analytic solutions to such a system however. Dec 16, 2015 at 18:15
• @george2079 , yes we are expecting almost 10 solutions. But not all of them are suitable for us. We already have found common solution for "$equ4$ and $equ5$""equ5 and equ6". So if we can find a common solution for "equ4 and equ6" then somehow by combining these invariants we will be able to find the solution for whole system. Dec 16, 2015 at 18:23
• on further though, you can trivially find infinitely many solutions. Just assume f a linear combination of the coefficients ( f=c1 x1 + c2 x2 .. ) You can arbitrarily specify three of the coefficients and algebraically solve. (If that's not satisfactory it indicates you have additional requirements you need to specify mathematically ) Dec 16, 2015 at 18:44
• If this set of equations has a solution, it seems likely that it will be an arbitrary function of three functions (probably polynomials) of the six variables, v1=x2 /x1, v2=x3 /x1, w1=y2/y1, w2=y3/y1, k1=z2/z1, k2=z3/z1. The challenge, then, is to find those three function. Dec 17, 2015 at 23:41

As noted in the question, this computation can be simplified by the substitution,

f[x1, x2, x3, y1, y2, y3, z1, z2, z3] := g[x2/x1, x3/x1, y2/y1, y3/y1, z2/z1, z3/z1]


in which case the six equations become

Simplify[{equ1, equ2, equ3}]
(* {0, 0, 0} *)

equ4 = Simplify[Simplify[equ4]/x1 /. {x2 -> x1 v2, x3 -> x1 v3, y2 -> y1 w2,
y3 -> y1 w3, z2 -> z1 k2, z3 -> z1 k3}]
(* (v2 + v3)*w3*Derivative[0, 0, 0, 1, 0, 0][g][v2, v3, w2, w3, k2, k3] +
v2*w2*Derivative[0, 0, 1, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] -
v3*(1 + 2*v2 + v3)*Derivative[0, 1, 0, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] -
v2*(1 + v2)*Derivative[1, 0, 0, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] *)

equ5 = Simplify[Simplify[equ5]/y1 /. {x2 -> x1 v2, x3 -> x1 v3, y2 -> y1 w2,
y3 -> y1 w3, z2 -> z1 k2, z3 -> z1 k3}]
(* k3*(w2 + w3)*Derivative[0, 0, 0, 0, 0, 1][g][v2, v3, w2, w3, k2, k3] +
k2*w2*Derivative[0, 0, 0, 0, 1, 0][g][v2, v3, w2, w3, k2, k3] -
w3*(1 + 2*w2 + w3)*Derivative[0, 0, 0, 1, 0, 0][g][v2, v3, w2, w3, k2, k3] -
w2*(1 + w2)*Derivative[0, 0, 1, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] +
v3*(1 + w2)*Derivative[0, 1, 0, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] +
v2*Derivative[1, 0, 0, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] *)

equ6 = Simplify[Simplify[ equ6]/z1 /. {x2 -> x1 v2, x3 -> x1 v3, y2 -> y1 w2,
y3 -> y1 w3, z2 -> z1 k2, z3 -> z1 k3}]
(* -(k3*(1 + 2*k2 + k3)*Derivative[0, 0, 0, 0, 0, 1][g][v2, v3, w2, w3, k2, k3]) -
k2*(1 + k2)*Derivative[0, 0, 0, 0, 1, 0][g][v2, v3, w2, w3, k2, k3] +
(1 + k2)*w3*Derivative[0, 0, 0, 1, 0, 0][g][v2, v3, w2, w3, k2, k3] +
w2*Derivative[0, 0, 1, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] *)


This system of first-order PDEs can be solved using the procedure described in Chapter V, Sec IV of Goursat's Differential Equations. The first step is to find the complete, non-commutative group of differential operators that includes equ4, equ5, and equ6. To do so, use the function comm, generalized from the related but simpler question, 99523.

drv = {Derivative[1, 0, 0, 0, 0, 0][g][v2, v3, w2, w3, k2, k3],
Derivative[0, 1, 0, 0, 0, 0][g][v2, v3, w2, w3, k2, k3],
Derivative[0, 0, 1, 0, 0, 0][g][v2, v3, w2, w3, k2, k3],
Derivative[0, 0, 0, 1, 0, 0][g][v2, v3, w2, w3, k2, k3],
Derivative[0, 0, 0, 0, 1, 0][g][v2, v3, w2, w3, k2, k3],
Derivative[0, 0, 0, 0, 0, 1][g][v2, v3, w2, w3, k2, k3]};
comm[equa_, equb_] := Collect[
(equa /. {Derivative[1, 0, 0, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] -> D[equb, v2],
Derivative[0, 1, 0, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] -> D[equb, v3],
Derivative[0, 0, 1, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] -> D[equb, w2],
Derivative[0, 0, 0, 1, 0, 0][g][v2, v3, w2, w3, k2, k3] -> D[equb, w3],
Derivative[0, 0, 0, 0, 1, 0][g][v2, v3, w2, w3, k2, k3] -> D[equb, k2],
Derivative[0, 0, 0, 0, 0, 1][g][v2, v3, w2, w3, k2, k3] -> D[equb, k3]}) -
(equb /. {Derivative[1, 0, 0, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] -> D[equa, v2],
Derivative[0, 1, 0, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] -> D[equa, v3],
Derivative[0, 0, 1, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] -> D[equa, w2],
Derivative[0, 0, 0, 1, 0, 0][g][v2, v3, w2, w3, k2, k3] -> D[equa, w3],
Derivative[0, 0, 0, 0, 1, 0][g][v2, v3, w2, w3, k2, k3] -> D[equa, k2],
Derivative[0, 0, 0, 0, 0, 1][g][v2, v3, w2, w3, k2, k3] -> D[equa, k3]}),
drv, Simplify]

equ7 = comm[equ4, equ5]
(* (k3*v2*w2 + k3*(v2 + v3)*w3)*Derivative[0, 0, 0, 0, 0, 1][g][v2, v3, w2, w3, k2, k3] +
k2*v2*w2*Derivative[0, 0, 0, 0, 1, 0][g][v2, v3, w2, w3, k2, k3] -
w3*(v3*(1 + w2 + w3) + v2*(1 + 2*w2 + w3))*
Derivative[0, 0, 0, 1, 0, 0][g][v2, v3, w2, w3, k2, k3] -
v2*w2*(1 + w2)*Derivative[0, 0, 1, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] +
v3*(v3*(1 + w2) + v2*(2 + w2))*Derivative[0, 1, 0, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] +
v2^2*Derivative[1, 0, 0, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] *)

equ8 = comm[equ5, equ6]
(* -(k3*((1 + 2*k2 + k3)*w2 + (1 + k2 + k3)*w3)*
Derivative[0, 0, 0, 0, 0, 1][g][v2, v3, w2, w3, k2, k3]) -
k2*(1 + k2)*w2*Derivative[0, 0, 0, 0, 1, 0][g][v2, v3, w2, w3, k2, k3] +
w3*((2 + k2)*w2 + (1 + k2)*w3)*Derivative[0, 0, 0, 1, 0, 0][g][v2, v3, w2, w3, k2, k3] +
w2^2*Derivative[0, 0, 1, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] -
v3*w2*Derivative[0, 1, 0, 0, 0, 0][g][v2, v3, w2, w3, k2, k3] *)


which are independent of the first three operators, increasing the size of the group to five. comm[equ4, equ6] vanishes identically and so does not add an operator. On the other hand, the seven additional commutators involving equ7 and equ8 yield expressions that are linear combinations of {equ4, equ5, equ6, equ7, equ8}. Thus, these five operators comprise the entire group.

From this information alone, we know that g is an arbitrary function of precisely one first integral. This first integral can be obtained by systematically eliminating variables and equations, one pair at a time, until a single equation of two variable remains. Start by solving any one of the equations. Because this process is presented in detail in 99523, here we merely list the required code without comment.

DSolve[equ4 == 0, g[v2, v3, w2, w3, k2, k3], {v2, v3, w2, w3, k2, k3}];
g[v2_, v3_, w2_, w3_, k2_, k3_] :=
h[w2, (v2 (1 + v2 + v3))/v3, (1 + v2) w2, (v3 w3)/v2, k2, k3];
tr1 = {equ4, equ5, equ6, equ7, equ8} // Simplify;
solw2 = First@tr1/(v2 w2) // Simplify;
newvar = Solve[Thread[{b1, b2, b3, b4, b5, b6} == List @@ solw2],
{v2, v3, w2, w3, k2, k3}] // Flatten;
tr1p = Collect[FullSimplify[Rest[tr1] /. solw2 -> 0 /. newvar], b1, FullSimplify]
/. b1*(z__) -> 0 /. Derivative[0, n1_, n2_, n3_, n4_, n5_][h][b1, b2, b3, b4, b5, b6]
-> Derivative[n1, n2, n3, n4, n5][h][b2, b3, b4, b5, b6];

DSolve[First@tr1p == 0, h[b2, b3, b4, b5, b6], {b2, b3, b4, b5, b6}] /. Log[z_] -> z;
h[b2_, b3_, b4_, b5_, b6_] := j[b4, b3, b5, (1 + b2 + b3) b4 b6, b6 (1 - b4)];
tr2 = tr1p // Simplify;
solb4 = First@tr2/(b4 (b4 - 1)) // Simplify;
newvar = Solve[Thread[{c1, c2, c3, c4, c5} == List @@ solb4], {b2, b3, b4, b5, b6}]
// Flatten;
tr2p = Collect[(Cancel[(c1 - 1) Rest@tr2] /. solb4 -> 0 /. newvar) //
FullSimplify, c1, FullSimplify];
tr2p[[3]] = tr2p[[3]]/c1;
tr2p = tr2p /. c1 z__ -> 0 /. Derivative[0, n1_, n2_, n3_, n4_][j][c1, c2, c3, c4, c5]
-> Derivative[n1, n2, n3, n4][j][c2, c3, c4, c5];

DSolve[Last@tr2p == 0, j[c2, c3, c4, c5], {c2, c3, c4, c5}];
j[c2_, c3_, c4_, c5_] := l[c5, c2, c3, (c2 - c4)/(1 + c3 + c5)];
tr3 = -tr2p // Simplify // RotateRight;
solc5 = First@tr3/(c5 (1 + c3 + c5)) // Simplify;
newvar = Solve[Thread[{d1, d2, d3, d4} == List @@ solc5], {c2, c3, c4, c5}] // Flatten;
tr3p = Collect[(Rest@tr3 /. solc5 -> 0 /. newvar) // FullSimplify, d1,
FullSimplify] /. d1 z__ -> 0 /. Derivative[0, n1_, n2_, n3_][l][d1, d2, d3, d4]
-> Derivative[n1, n2, n3][l][d2, d3, d4];

DSolve[Last@tr3p == 0, l[d2, d3, d4], {d2, d3, d4}] // Simplify;
l[d2_, d3_, d4_] := m[d3, (1 + d2) d3, (d2 (1 + d2 - d4))/d4];
tr4 = tr3p // Simplify // RotateRight;
sold3 = First@tr4/(d2 d3);
newvar = Solve[Thread[{e1, e2, e3} == List @@ sold3], {d2, d3, d4}] //  Flatten;
tr4p = Collect[(-(e2/e1) Rest@tr4 /. sold3 -> 0 /. newvar) //
FullSimplify, e1, FullSimplify] /. e1 z__ -> 0
/. Derivative[0, n1_, n2_][m][e1, e2, e3] -> Derivative[n1, n2][m][e2, e3];

(DSolve[Last@tr4p == 0, m[e2, e3], {e2, e3}] // Simplify) /. Log[z_] -> z;
(((((((%[[1, 1, 2]] /. Thread[{e1, e2, e3} -> List @@ sold3]) // Simplify)
/. Thread[{d1, d2, d3, d4} -> List @@ solc5]) // Simplify)
/. Thread[{c1, c2, c3, c4, c5} -> List @@ solb4]) // Simplify)
/. Thread[{b1, b2, b3, b4, b5, b6} -> List @@ solw2] // Simplify)
/. {v2 -> x2/x1, v3 -> x3/x1, w2 -> y2/y1, w3 -> y3/y1, k2 -> z2/z1, k3 -> z3/z1})
// Simplify

(* C[1][(((x2 y2 z2 + x1 (y2 z2 + y1 (z1 + z2))) (x3 y3 z3 + x2 (y3 z3 + y2 (z2 + z3))))/
(x2 y2 z2 (x3 y3 z3 + x2 (y3 z3 + y2 (z2 + z3)) +
x1 (y3 z3 + y2 (z2 + z3) + y1 (z1 + z2 + z3)))))] *)

• Oh! my friend. You found the solution!! I checked it! it satisfied! I am very very happy :) Thank you very much! :) Can you please let me know how did you found it? Dec 26, 2015 at 12:18
• @farrokh I shall provide the derivation within a few hours. For now, see your earlier question, which I also solved. Dec 26, 2015 at 12:19
• Oh! I am very happy :) Thanks a lot! :) Dec 26, 2015 at 12:19