# Numerically integrating exact differentials

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 + 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 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?

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

• Thank you. My actual PDE is a bit more complicated than the one I gave in the example, so I need a numerical approach. Your numerical solution is very nice, but it does involve a symbolic Integrate step which won't work for me. If that could be changed to an NIntegrate step somehow, then that would be very nice. – Bruce Bartlett Feb 13 at 8:52
• @BruceBartlett NIntegrate computes only definite integrals. It's better to use NDSolve to compute the antiderivative. Please see my update. – Michael E2 Feb 13 at 14:35
• Ok, thank you, this works well. – Bruce Bartlett Feb 18 at 8:59

{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*)

• Thank you, but the actual problem that I have does not admit an analytic solution by DSolve, so it can't be solved step by step in this way. I just used this simple example as a base case. – Bruce Bartlett Feb 13 at 8:49
• Did you read my answer? I did not use DSolve because your problem is not a PDE. You have the derivatives explicitly! – Daniel Huber Feb 13 at 8:53
• You are still using Integrate, i.e. a symbolic method. I need a purely numeric method. – Bruce Bartlett Feb 13 at 12:22
• No problem, simply use NIntegrate. – Daniel Huber Feb 13 at 13:01
• Ok thank you. Your answer has many useful features. – Bruce Bartlett Feb 18 at 8:59