Numerically integrating exact differentials

Question

I would like to be able to numerically integrate an exact differential 1-form. Mathematica can do this symbolically using DSolve (though I don't know how to insert an initial condition), but I get a "system is overdetermined" error when I try to do this numerically using NDSolve.

For instance, consider the following symbolic example posted by bbgodfrey on this forum:

DSolveValue[
   {D[f[x, y], x] == 1 + 2 E^(x^2 + y^2) x Cos[x + y] - E^(x^2 + y^2) Sin[x + y], 
    D[f[x, y], y] ==  3 y^2 + 2 E^(x^2 + y^2) y Cos[x + y] - E^(x^2 + y^2) Sin[x + y]}, 
    f[x, y], {x, y}] // FullSimplify

output: x + y^3 + C[1] + E^(x^2 + y^2) Cos[x + y]

That works fine... but there are two problems. The first (smaller) problem is that putting in the initial condition f(0,0)=0 (which forces c[1] to be -1) doesn't work:

DSolveValue[{D[f[x, y], x] == 
    1 + 2 E^(x^2 + y^2) x Cos[x + y] - E^(x^2 + y^2) Sin[x + y], 
   D[f[x, y], y] == 
    3 y^2 + 2 E^(x^2 + y^2) y Cos[x + y] - E^(x^2 + y^2) Sin[x + y], 
   f[0, 0] == 0}, f[x, y], {x, y}] // FullSimplify

output: DSolveValue[{D[f[x, y], x] == 
    1 + 2 E^(x^2 + y^2) x Cos[x + y] - E^(x^2 + y^2) Sin[x + y], 
   D[f[x, y], y] == 
    3 y^2 + 2 E^(x^2 + y^2) y Cos[x + y] - E^(x^2 + y^2) Sin[x + y], 
   f[0, 0] == 0}, f[x, y], {x, y}] // FullSimplify

But the second (more important) problem is that I can't get NDSolve to numerically solve the exact PDE (with the initial condition f(0,0)=0):

NDSolveValue[{D[f[x, y], x] == 
   1 + 2 E^(x^2 + y^2) x Cos[x + y] - E^(x^2 + y^2) Sin[x + y], 
  D[f[x, y], y] == 
   3 y^2 + 2 E^(x^2 + y^2) y Cos[x + y] - E^(x^2 + y^2) Sin[x + y], 
  f[0, 0] == 0}, f, {x, -2, 2}, {y, -2, 2}]

output: NDSolveValue::overdet: There are fewer dependent variables, {f[x,y]}, than equations, so the system is overdetermined.

Can anyone help?

Answer 1

The problem is that an exact system is an overdetermined yet consistent system and so has a special solution. This can be seen in the initial condition: Normally it is defined along a line f[0, y] == ..., not just at a point. We don't have enough degrees of freedom to solve a condition along a line.

Recently (V12.2) the problem domain of DSolve was extended to include exact PDEs. I am not sure why DSolve does not solve them with an initial point-condition, but apparently it does not. It may be they are still working on reliably parsing the problem. Consequently we seem to have to solve for the constant with separate code.

The NDSolve problem can be alleviated if one of the PDEs can be integrated and used as an initial condition. (Updated to numerically integrate the initial condition.)

dsol = First@
   DSolve[{D[f[x, y], x] == 
      1 + 2 E^(x^2 + y^2) x Cos[x + y] - E^(x^2 + y^2) Sin[x + y], 
     D[f[x, y], y] == 
      3 y^2 + 2 E^(x^2 + y^2) y Cos[x + y] - E^(x^2 + y^2) Sin[x + y]},
    f, {x, y}];
csol = First@Solve[f[0, 0] == 0 /. dsol, C[1]];
fSOL = f[x, y] /. dsol /. csol; (* particular solution *)

fIFN = NDSolveValue[
  {D[f[x, y], x] == 
    1 + 2 E^(x^2 + y^2) x Cos[x + y] - E^(x^2 + y^2) Sin[x + y]
   , f[0, y] == # - (# /. y -> 0) &@(*to enforce IC*)
       NDSolveValue[
         {F'[y] == 3 y^2 + 2 E^(x^2 + y^2) y Cos[x + y] - 
             E^(x^2 + y^2) Sin[x + y] /. x -> 0, F[0] == 0}
         , F[y]
         , {y, -2, 2}
         ]},
  f,
  {x, -2, 2}, {y, -2, 2}]

Plot3D[{fSOL, fIFN[x, y]}, {x, -2, 2}, {y, -2, 2}]

The blotchy color indicates the surfaces nearly coincide.

Appendix

Original numerical method symbolically integrated the boundary condition:

fIFN = NDSolveValue[{D[f[x, y], x] == 
     1 + 2 E^(x^2 + y^2) x Cos[x + y] - E^(x^2 + y^2) Sin[x + y],
    f[x, y] == # - (# /. y -> 0) &@ (* to enforce IC *)
      Integrate[
       3 y^2 + 2 E^(x^2 + y^2) y Cos[x + y] - 
        E^(x^2 + y^2) Sin[x + y], y] /. x -> 0}, 
   f, {x, -2, 2}, {y, -2, 2}];

Answer 2

Consider your equations:

{D[f[x, y], x] == 1 + 2 E^(x^2 + y^2) x Cos[x + y] - E^(x^2 + y^2) Sin[x + y],

D[f[x, y], y] ==  3 y^2 + 2 E^(x^2 + y^2) y Cos[x + y] - E^(x^2 + y^2) Sin[x + y]}

That's not actually a PDE!. You have the derivatives as explicit functions. We can get the function by simple integration, where you have to take care of the addition constant (or better addition function here).Towards this aim, we integrate D[f[x, y], x] and D[f[x, y], y], call them fdx and fdy and get s1 and s2:

fdx[x_] = 1 + 2 E^(x^2 + y^2) x Cos[x + y] - E^(x^2 + y^2) Sin[x + y];
s1 = Integrate[fdx[x], x];

fdy[y_] = 3 y^2 + 2 E^(x^2 + y^2) y Cos[x + y] - E^(x^2 + y^2) Sin[x + y];
s2 = Integrate[fdy[y], y]

Taking care of the constant of integration, we must add to s1 a function of only y, call it fy. And to s2 a function of only x, call it fx. Then we should get the same function, the difference should be zero:

(s1 + fy) - (s2 + fx) == 0 // Simplify
(* fx + y^3 == fy + x *)

or

fx == x + const  and fy== y^3 + const

Therefore the searched for function f is:

f == s1+ y^3 + const == s2 + x + const

with some constan const, determine by the initial condition. We can check if both definitions are the same:

s1 + y^3 + const == s2 + x + const // Simplify
(*True*)

To determine const we use e.g. the initial condition f[0,0]==0;

s1 + y^3 + const == 0 /. {x -> 0, y -> 0}
(* 1 + const == 0 *)

Therefore, const== -1 and finally f:

f[x_, y_] = s1 + y^3 - 1

what is:

Finally we can check if this function really fulfills the original equations:

{D[f[x, y], x] == 1 + 2 E^(x^2 + y^2) x Cos[x + y] - E^(x^2 + y^2) Sin[x + y], 
  D[f[x, y], y] == 3 y^2 + 2 E^(x^2 + y^2) y Cos[x + y] -  E^(x^2 + y^2) Sin[x + y]} // Simplify

(* True, True*)