Solving an 'odd' differential equation with NDSolve

I need to solve a differential equation of the type

$\qquad \partial_{x_1}y(x_1,x_2)= y(x_1,x_1)\,y(x_1,x_2)$

with initial condition

$\qquad y(0,x_2)=x_2$.

Now if I try to solve this with NDSolve I get an error that tells me "the arguments should be ordered consistently" (NDSolve::conarg), here is my code:

s =
NDSolve[
{D[y[x1, x2], x1] == y[x1, x1] y[x1, x2], y[0, x2] == x2},
y, {x1, 0, 1}, {x2, 0, 1}]


On my daylong search I haven't found any helpful solution for this seemingly easy problem.

• Since it's not mentioned in the question, I'll point out that the cause of the error is that the arguments of y[x1, x1] and y[x1, x2] are not the same. Sep 10 '18 at 16:35
• The y[x1,x1] is the point of the problem, otherwise it would be trivial. Sorry, I should have clarified that. I want it to work this way. Sep 10 '18 at 17:16
• Just wanted to make sure the meaning of the error was clear. I wonder, what makes this problem seem easy? Is it a standard type? Are there standard methods for numerically integrating it? Sep 10 '18 at 22:09
• Welcome to Mathematica.SE! I hope you will become a regular contributor. 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. Sep 12 '18 at 3:55

This PDE can be solved symbolically. To begin, replace y[x1, x1] by f[x1], which DSolve can handle without difficulty.

DSolve[{D[y[x1, x2], x1] == f[x1] y[x1, x2], y[0, x2] == x2},
y[x1, x2], {x1, x2}] [[1, 1]]
(* y[x1, x2] -> E^(Integrate[-f[K[1]], {K[1], 1, 0}]
- Integrate[-f[K[1]], {K[1], 1, x1}])*x2


From this solution, it is clear that y[x1, x2]/x2 is independent of x2. Consequently, y[x1, x1]/x1 == y[x1, x2]/x2. Substituting this identity into the PDE then yields,

s = DSolve[{D[y[x1, x2], x1] == (x1/x2) y[x1, x2]^2, y[0, x2] == x2},
y[x1, x2], {x1, x2}] [[1, 1]]
(* y[x1, x2] -> -((2 x2)/(-2 + x1^2)) *)

Plot3D[s//Last, {x1, 0, 1}, {x2, 0, 1}, AxesLabel -> {x1, x2, y},
ImageSize -> Large, LabelStyle -> {Bold, Black, Medium}]


Not surprisingly, the plot closely resembles the numerical result of Michael E2. Note that this method works with this PDE for any initial condition.

This solution can, if desired, be validated by back-substitution:

{D[y[x1, x2], x1] == y[x1, x1] y[x1, x2], y[0, x2] == x2} /.
y -> Function[{x1, x2}, -((2 x2)/(-2 + x1^2))]
(* {True, True} *)

• Writing down f[x1], you implicitly claim it a known function, but it is not so. Therefore, your solution is built on the sand. Am I not right? Sep 12 '18 at 4:57
• @user64494 No. I do not claim that f[x1] is known, only that y[x1, x1] is a function only of x1. Sep 12 '18 at 5:02
• OK, Function[{x1, x2}, -((2 x2)/(-2 + x1^2))] is a solution. I still have doubts concerning your DSolving. Sep 12 '18 at 5:20

One could try a iterative approach approximating y[x1, x1] with the previous, if the iterations converge to solution:

foo::div = "Error increased. Divergent?";
foo::slwcon = "Slow convergence. Error  after  steps.";
tol = 10^-6;
NestWhile[
Function[{s1, s2, err1, err2, iter},
With[{s =
NDSolveValue[{D[y[x1, x2], x1] == s1[x1, x1] y[x1, x2],
y[0, x2] == x2}, y, {x1, 0, 1}, {x2, 0, 1}]},
{s, s1,
Flatten[s1["ValuesOnGrid"], 1] - s @@@ Flatten[s1["Grid"], 1] //
Abs // Max,
err1,
iter + 1}
]] @@ # &,
{NDSolveValue[{D[y[x1, x2], x1] == 1*y[x1, x2], y[0, x2] == x2},
y, {x1, 0, 1}, {x2, 0, 1}],
1 &,
10., 100., 0},
Function[{s1, s2, err1, err2, iter},
If[err1 > err2,
Message[foo::div]; False,
If[err1 > tol,
If[iter == 10, Message[foo::slwcon, err1, iter]];
True,
False]
]
] @@ # &,
1, 100
]
s = {y -> First[%]};


Plot3D[y[x1, x2] /. s, {x1, 0, 1}, {x2, 0, 1}, AxesLabel -> Automatic]


The error is so-so:

Plot3D[D[y[x1, x2], x1] - y[x1, x1] y[x1, x2] /. s // Evaluate,
{x1, 0, 1}, {x2, 0, 1}, AxesLabel -> Automatic]


• Thanks for your effort! Sep 11 '18 at 14:38