# 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
Commented Oct 28, 2018 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? Commented Oct 28, 2018 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. Commented Oct 28, 2018 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. Commented Oct 28, 2018 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. Commented Oct 29, 2018 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. Commented Oct 29, 2018 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. Commented Oct 29, 2018 at 13:24
• You are right. Thank you for the insight though. Commented Oct 29, 2018 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. Commented Oct 29, 2018 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? Commented Oct 29, 2018 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. Commented Oct 29, 2018 at 13:02
• But this is a completely different task. In which area is it necessary to minimize the functional R — finite or infinite? Commented Oct 29, 2018 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. Commented Oct 29, 2018 at 13:27