# How can I solve a system of PDEs with constraints but having unknown boundary conditions

I have a 1st order system of two PDEs with two independent variables and two dependent variables.

mu = -0.1;
lambda = -1;
x10 = 0;
x20 = 0;
eq =
{mu*x1 + D[h1[x1, x2], x1]*mu*x1 + D[h1[x1, x2], x2]*lambda*(x2 - x1^2) ==
mu*(x1 + h1[x1, x2]),
lambda*(x2 - x1^2) + D[h2[x1, x2], x1]*mu*x1 + D[h2[x1, x2], x2]*lambda*(x2 - x1^2) ==
lambda*(x2 + h2[x1, x2])};


The domain is given by {x1, -2, 2} and {x2, -2, 2}. The problem is that I don't know the boundary condition of the PDE, but I do have a constraint that the gradient of h1 and h2 at point {x10, x20} is zero.

D[h1[x1, x2], x1]/.{x1 -> x10, x2 -> x20} == 0
D[h1[x1, x2], x2]/.{x1 -> x10, x2 -> x20} == 0
D[h2[x1, x2], x1]/.{x1 -> x10, x2 -> x20} == 0
D[h2[x1, x2], x2]/.{x1 -> x10, x2 -> x20} == 0


How can I compute an approximation to h1[x1, x2] and h2[x1, x2] numerically?

• That is not enough information. For a numerical integration you need a definite region. And you need initial and boundary conditions. Derivatives at only one point are not enough. Commented Jan 17, 2021 at 8:58
• @DanielHuber The region is {x1, -2, 2} and {x2, -2, 2}. But the initial and boundary conditions are unknown. Commented Jan 17, 2021 at 9:09
• Then there is no unique solution. Commented Jan 17, 2021 at 10:33
• This needs work. First this is not a system of coupled PDE, you can solve it separately. Further, look at your equations for x1==x2==0. The equation for h1 then reads: 0. == -0.1 h1[0, 0] what does not define the derivatives. Likewise for h2. Commented Jan 17, 2021 at 13:32
• I’m voting to close this question because the issue it raises is not really a Mathematica issue but a matter of the OP not having grasped the relevant mathematics. Commented Jan 17, 2021 at 19:00

The two PDE are independent, and solutions to the second can be obtained without difficulty. First, for convenience, move all terms to the left side of the equation and also rationalize mu.

Subtract @@ eq[[2]]
(* x1^2 + h2[x1, x2] - (-x1^2 + x2)*D[h2[x1, x2],x2] - (x1*D[h2[x1, x2],x1])/10 *)


x2 appears explicitly only as the coefficient of D[h2[x1, x2],x2], so let us suppose that h2 is a function of x1 only. (This may exclude some solutions, but the OP asked for any solution.)

% /. h2 -> Function[{x1, x2}, h2[x1]]
(* x1^2 + h2[x1] - 1/10 x1 h2'[x1] *)

DSolveValue[% == 0, h2[x1], x1]
(* -((5 x1^2)/4) + x1^10 C[1] *)


It also satisfies D[h2[x1, x2], x1]/.{x1 -> x10, x2 -> x20} == 0,as desired.

A similar calculation can be performed for the first equation.

Subtract @@ eq[[1]];
% /. h1 -> Function[{x1, x2}, h1[x1]];
DSolveValue[% == 0, h1[x1], x1]
(* x1 C[1] *)


However, the requirement that D[h1[x1, x2], x1]/.{x1 -> x10, x2 -> x20} == 0 requires that C[1] = 0, yielding only the trivial solution.

Further investigation shows that DSolve can provide general solutions for the two PDEs.

s1 = DSolveValue[eq[[1]], h1[x1, x2], {x1, x2}]
(* x1 C[1][-(5/(4 x1^8)) + x2/x1^10] *)

s2 = DSolveValue[eq[[2]], h2[x1, x2], {x1, x2}]
(* -((5 x1^2)/4) + x1^10 C[1][-(5/(4 x1^8)) + x2/x1^10] *)


where C[1] is an arbitrary function of z = (x2 - 5 x1^2/4)/x1^10. (In fact, DSolve cans solve many quasilinear PDEs.) To apply the constraints at {x10, x20}, C[1] must be well behaved there, as it is for the specific solutions derived earlier in this answer (for which C[1] is constant). Possibly, other specific solutions satisfying the constraints can be found among the many possible C[1] functions. For instance Sinh[z]/Cosh[z]^3 is continuously differentiable at (0, 0} and satisfies the constraints, although it is not analytical there. Whether such a solution is acceptable depends of how the OP plans to use it.

The OP posed a similar but coupled pair of PDEs in a comment below:

eq1 = {-1/10*x1 - (-x1^2 + x2)*D[h1[x1, x2], x2] - x1*D[h1[x1, x2], x1]/10 ==
(-x1 - h1[x1, x2])/10 - h2[x1, x2],
x1^2 - x2 - (-x1^2 + x2)*D[h2[x1, x2], x2] - x1*D[h2[x1, x2], x1]/10 ==
-x2 - h2[x1, x2]}


DSolve can solve these ODEs, but only with assistance. Because eq1[[2]] is as before, its solution determined above can be inserted into eq1[[1]].

eq1s = Simplify[Subtract @@ eq1[[1]] /.
h2 -> Function[{x1, x2}, x1^10 c1[(x2 - 5 x1^2/4)/x1^10] - 5 x1^2/4]]
(* (-5*x1^2)/4 + x1^10*c1[((-5*x1^2)/4 + x2)/x1^10] + h1[x1, x2]/10 +
x1^2*D[h1[x1, x2], x2] - x2*D[h1[x1, x2], x2] - x1*D[h1][x1, x2], x1]/10 *)


The inhomogeneous term -((5 x1^2)/4) contributes -(25 x1^2)/2 to the solution for h1, as can be seen from evaluating

DSolve[-((5*x1^2)/4) + (1/10)*h1[x1, x2] - (1/10)*x1*D[h1[x1, x2], x1] == 0,
h1[x1, x2], x1]


On this basis perform the second substitution,

Simplify[eq1s /.
h1 -> Function[{x1, x2}, h1[x1] c c1[(x2 - 5 x1^2/4)/x1^10] - 25 x1^2/2]]
(* 1/10 c1[(-((5 x1^2)/4) + x2)/x1^10] (10 x1^10 + c h1[x1] - c x1 D[h1[x1], x1]) *)


(h1 and h2 must involve the same function c1, up to a constant multiplier designated c, in order to solve the two equations. The solution to this last ODE is

DSolveValue[% == 0, h1[x1], x1]
(* (10 x1^10)/(9 c) + x1 C[1] *)


Gathering the solutions together and back-substituting them into eq1 verifies their validity.

Simplify[eq1 /.
{h1 -> Function[{x1, x2}, (10 x1^10/9 + c x1) c2[(x2 - 5 x1^2/4)/x1^10] - 25 x1^2/2],
h2 -> Function[{x1, x2}, x1^10 c2[(x2 - 5 x1^2/4)/x1^10] - 5 x1^2/4]}]
(* {True, True} *)


For the case of c2 constant and the constraints applied, the solution becomes

h1 -> Function[{x1, x2}, 10 x1^10/9 c2 - 25 x1^2/2]
h2 -> Function[{x1, x2}, x1^10 c2 - 5 x1^2/4]}]

• Hi, thanks for your answer. In general, the two PDEs are not independent. Just this simple example decouples the two PDEs. The solution given in the Addendum is very interesting. Will it work if the two PDEs are dependent on each other? For example, the PDEs look like this: eq = {mu*x1 + D[h1[x1, x2], x1]*mu*x1 + D[h1[x1, x2], x2]*lambda*(x2 - x1^2) == mu*(x1 + h1[x1, x2])+lambda*h2[x1,x2], lambda*(x2 - x1^2) + D[h2[x1, x2], x1]*mu*x1 + D[h2[x1, x2], x2]*lambda*(x2 - x1^2) == lambda*(x2 + h2[x1, x2])}; Commented Jan 18, 2021 at 7:13
• @LaZheng In general, coupling the two PDEs makes them much more difficult to solve. For the equations in your comment, however, it should not be difficult. I shall do so later today. Commented Jan 18, 2021 at 16:16
• @LeZHENG As you requested I have added the general solution to the coupled equations in your comment. This was relatively easy, because the two PDEs still had the same Characteristic surfaces even after coupling. Note that the general solution is not analytic at {0, 0}, which may be a consideration in choosing which solutions to use. The exception of course is to take c2 as a constant. By the way, I voted to reopen your question, but that probably won't happen. Too bad. Commented Jan 19, 2021 at 23:56
• @LeZHENG I hope you will become a regular contributor to Mathematica.SE. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Best wishes. Commented Jan 20, 2021 at 0:05
• Never mind, maybe it's really bad question and your answer is a new year gift for me. It's the first time that I see what DSolveValue can do. And I got another PDE to solve, the format is like the above, but maybe because it's complicated that DSolveValue gives no answer. I post the codes in the comment below. Commented Jan 20, 2021 at 1:03