Newbie PDE 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.

• NDSolve won't solve for general solutions. Add exactly enough initial conditions for your system to specify a specific solution and see what happens. – Bill Oct 28 '18 at 17:13
• @Bill well I tried to add a lot of symmetry and initial conditions. It seems never enough. How many do I need? – Cowboy Trader Oct 28 '18 at 17:34
• Apparently, we need to develop an algorithm for solving nonlinear equations in 3D. In this case, we have an obvious message: NDSolve::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve. – Alex Trounev Oct 28 '18 at 21:07
• If there is an interest in a numerical solution in 3D, I can suggest an algorithm that I use to solve the Navier-Stokes equation. – Alex Trounev Oct 28 '18 at 22:15

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

• It is for certain that solution of c1 is not only a function of y1. However there are symmetries such as c1[y1,y2,y3] == c1[y1,y3,y2], c2[y1,y2,y3] ==c2[y3,y2,y1], c3[y1,y2,y3]==c3[y2,y1,y3], which holds the position of yn in cn and changes the other two (I think these symmetries will apply to the derivatives as well). I am trying to understand what you did. "/. C[1] -> c" don't know what this expression means, new to Mathematica. – Cowboy Trader Oct 29 '18 at 7:13
• @CowboyTrader /. C[1] -> c renames the integration constant C[1] as the constant c. This is merely for notational simplicity and has no substantive impact on the solution. By the way, the solution I derived satisfies the symmetries listed in your comment. – bbgodfrey Oct 29 '18 at 13:13
• @CowboyTrader In response to your comment to Alex Trounev, I do not believe that it is possible to minimize R everywhere. Setting Grad[R] to zero instead specifies that R is a constant. So, I suggest you focus on correctly formulating your mathematical problem, and only then thinking about how to solve it with Mathematica or some other means. Best wishes. – bbgodfrey Oct 29 '18 at 13:24
• You are right. Thank you for the insight though. – Cowboy Trader Oct 29 '18 at 13:30

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


• It does not matter which constant to use in the initial data. You have a correctly posed problem and a method of numerical solution. Change the initial data and get the solution that you need. – Alex Trounev Oct 29 '18 at 11:56
• The numerical model will not change with variable R=R[y1,y2,y3]. However, the requirement of a zero gradient in any sense means that R is constant. Explain the origin of the equation. Is this R a curvature? – Alex Trounev Oct 29 '18 at 12:54
• Basically I am trying to find the functions c1,c2,c3 that minimize R for every possible point. y1,y2,y3. I don't see why this should mean R should be a constant. – Cowboy Trader Oct 29 '18 at 13:02
• But this is a completely different task. In which area is it necessary to minimize the functional R — finite or infinite? – Alex Trounev Oct 29 '18 at 13:22
• -infinity to +infinity. Sorry to waste your precious time. I am just new to this. I tried functional derivative as well (with VariationalD function). But it also leads to a constant R for some reason. – Cowboy Trader Oct 29 '18 at 13:27