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?
ClearAll["Global`*"]
<< 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.