Newbie PDE Question

Question

I downloaded Mathematica trial in order to try to solve a PDE that Matlab symbolic couldn't handle. I don't know all the bells and whistles but going along with online documentation.

I am trying to solve a set of differential equations. I typed the following:

R = (y1 - c1[y1, y2, y3])^2 + (y2 - c2[y1, y2, y3])^2 + (y3 - c3[y1, y2, y3])^2 + 
  2*(D[c1[y1, y2, y3], y1] + D[c2[y1, y2, y3], y2] + D[c3[y1, y2, y3], y3])

DSolve[{D[R, y1] == 0, D[R, y2] == 0, D[R, y3] == 0}, {c1[y1, y2, y3],
   c2[y1, y2, y3], c3[y1, y2, y3]}, {y1, y2, y3}]

DSolve is giving my equation back to me without any error messages. I thought maybe there is no analytical solution so I tried NDSolve, as well. NDSolve complains that the system is underdetermined. I threw in a lot of symmetry conditions and an initial value with no avail.

These are the symmetry conditions, which I think apply to derivatives as well.

c1[y1,y2,y3] == c1[y1,y3,y2], c2[y1,y2,y3] ==c2[y3,y2,y1], c3[y1,y2,y3]==c3[y2,y1,y3]

As a newbie I must be doing something wrong, as it doesn't like my prose. Any pointers will be appreciated.

Answer 1

A solution can be obtained as follows. First, put aside for a moment the particular definition of R and, instead, write it as an arbitrary function of the three independent variables, R[y1, y2, y3]. DSolve can handle the resulting system without difficulty.

Clear[R]
Flatten@DSolve[{D[R[y1, y2, y3], y1] == 0, D[R[y1, y2, y3], y2] == 0, 
    D[R[y1, y2, y3], y3] == 0}, R[y1, y2, y3], {y1, y2, y3}] /. C[1] -> c
(* {R[y1, y2, y3] -> c} *)

where c is a constant. Not surprisingly, inserting this result into

DSolve[R == c, {c1[y1, y2, y3], c2[y1, y2, y3], c3[y1, y2, y3]}, {y1, y2, y3}]

fails with the error message

DSolve::underdet: There are more dependent variables than equations, so the system is underdetermined.

So, let us make the plausible assumption that c1 is a function of y1 only, and similarly for c2 and c3. Then, R == c can be decomposed into three ODEs, and c1[y1] is determined by

DSolve[cy1 == (y1 - c1[y1])^2 + 2 D[c1[y1], y1], c1[y1], y1] // Flatten
(* {c1[y1] -> Sqrt[-2 + cy1] + y1 + 
    2/(-(1/Sqrt[-2 + cy1]) - E^(Sqrt[-2 + cy1] y1) C[1])} *)

c2[y2] and c3[y3] are given by the same expression but with y1 and the constant cy1 replace by y2 and cy2, and y3 and cy3, respectively. Of course,

c == cy1 + cy2 + cy3

Answer 2

Equations {D[R, y1] == 0, D[R, y2] == 0, D[R, y3] == 0} mean that R is a constant. Therefore, one equation needs to be solved R == R0 = const. Consider a potential field {c1,c2,c3}=Grad[f[y1,y2,y3],{y1,y2,y3}], then the equation is reduced to the form $(\vec y-\nabla f)^2+2\nabla ^2f=R0$, $\vec y=(y1,y2,y3)$

This nonlinear equation in 3D cannot be solved even numerically with the help of NDSolve of Mathematica 11.3. Since Mathematica 11.3 implements the finite element method for solving linear problems for elliptic equations in 3D, we use this method together with the fixed point method for solving a nonlinear problem:

R0=6; F[0][y1_, y2_, y3_] := (y1^2 + y2^2 + y3^2)/2
Do[F[i] = NDSolveValue[{(y3 - 
\!\(\*SuperscriptBox[\(F[\(-1\) + i]\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[y1, y2, y3])^2 + (y2 - 
\!\(\*SuperscriptBox[\(F[\(-1\) + i]\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[y1, y2, y3])^2 + (y1 - 
\!\(\*SuperscriptBox[\(F[\(-1\) + i]\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[y1, y2, y3])^2 + 2 (
\!\(\*SuperscriptBox[\(f\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "2"}], ")"}],
Derivative],
MultilineFunction->None]\)[y1, y2, y3] + 
\!\(\*SuperscriptBox[\(f\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "2", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[y1, y2, y3] + 
\!\(\*SuperscriptBox[\(f\), 
TagBox[
RowBox[{"(", 
RowBox[{"2", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[y1, y2, y3]) == R0, 
     DirichletCondition[f[y1, y2, y3] == 0, True]}, 
    f, {y1, 0, 1}, {y2, 0, 1}, {y3, 0, 1}], {i, 1, k}];

The solution of the equation at the last step and the difference of solutions at the last two steps

{Plot3D[F[k][0.5, y2, y3], {y2, 0, 1}, {y3, 0, 1}, Mesh -> None, 
 ColorFunction -> Hue], Plot[F[k][.5, .5, y3] - F[k - 1][.5, .5, y3], {y3, 0, 1}]}

The distribution of the vector field {c1,c2,c3} in space, on the plane and the line level of one component c1[y1,y2,.5]

{VectorPlot3D[
  Evaluate[{D[F[k][y1, y2, y3], y1], D[F[k][y1, y2, y3], y2], 
    D[F[k][y1, y2, y3], y3]}], {y1, 0, 1}, {y2, 0, 1}, {y3, 0, 1}, 
  VectorColorFunction -> Hue], 
 VectorPlot[
  Evaluate[{D[F[k][y1, y2, y3], y1], D[F[k][y1, y2, y3], y2]} /. 
    y3 -> .5], {y1, 0, 1}, {y2, 0, 1}, VectorColorFunction -> Hue], 
 ContourPlot[
  Evaluate[D[F[k][y1, y2, y3], y1] /. y3 -> .5], {y1, 0, 1}, {y2, 0, 
   1}, Contours -> 20, ColorFunction -> Hue]}