In a Newton iterative FEM PDE solver I am writing, each iterative step involves updating the current solution function u by adding to it the increment function du, which is the solution of some linearised PDE whose coefficients depend on the current u. Both u and du are InterpolatingFunction.

Now, as the number of iterates n grows, it becomes quite inefficient to evaluate the sum of n interpolating functions.

How can two InterpolatingFunction over a finite element mesh be added together to get a single InterpolatingFunction, whilst keeping all the original unstructured mesh properties / interpolation order?

<< NDSolve`FEM`
rad = 5; th = Pi/4;
reg = Disk[{0,0}, rad, {-th,th}];
mesh = ToElementMesh[reg];

Say the two functions to be added are

f1 = NDSolveValue[{Laplacian[u[x,y],{x,y}] == 1^2 * u[x,y]
       + NeumannValue[1, True]}, u, Element[{x,y}, mesh]];
f2 = NDSolveValue[{Laplacian[u[x,y],{x,y}] == 2^2 * u[x,y]
       + NeumannValue[1, True]}, u, Element[{x,y}, mesh]];

Attempt 0.

sum0 = Function[{x,y}, f1[x,y] + f2[x,y] // Evaluate];
Do[  f1[0,0], 100000] // AbsoluteTiming  (* {0.58909, Null} *)
Do[sum0[0,0], 100000] // AbsoluteTiming  (* {1.64215, Null} *)

Using Function obviously doesn't lump the sum into a single InterpolatingFunction.

Attempt 1.

sum1 = NDSolveValue[s[x,y] == f1[x,y] + f2[x,y], s, Element[{x,y}, mesh]];

NDSolveValue::derivs: No derivatives of dependent variables were found in the equations. NDSolveValue is designed to solve differential or differential algebraic equations. Use NSolve or FindRoot to numerically solve algebraic equations.

See (23416); it's a pity you can't solve "zeroth order differential equations" anymore.

Attempt 2.

sum2 = FunctionInterpolation[f1[x,y] + f2[x,y], Element[{x,y}, mesh]];

FunctionInterpolation::range: "Argument {x,y}[Element]mesh is not in the form of a range specification, {x, xmin, xmax}."

Using the form {x, xmin, xmax}, ... is not an option because the reg is not rectangular and mesh is unstructured, and I want the sum to inherit these properties.


2 Answers 2


You can do something like this (assuming of course that the InterpolatingFunction f1 and f2 are generating using the same mesh):

f12 = ElementMeshInterpolation[{mesh}, f1["ValuesOnGrid"] + f2["ValuesOnGrid"]]

Where ElementMeshInterpolation is from the



f1[1, 1] + f2[1, 1]
(* -1.84159 *)

f12[1, 1]
(* -1.84159 *)

Do[f1[0, 0], 100000] // AbsoluteTiming
(* {0.223145, Null} *)

Do[f12[0, 0], 100000] // AbsoluteTiming
(* {0.213769, Null} *)

Attempt 3.

It turns out you can trick Mathematica into not complaining about lack of derivatives by adding a dependent derivative in a separate equation:

sum3 = NDSolveValue[
         {Laplacian[dummy[x, y], {x, y}] == 0,
          s[x, y] == f1[x, y] + f2[x, y]
         }, s, Element[{x, y}, mesh]]

Of course this is slower than chuy's method.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.