# How can I improve and accelerate my 2 D algorithm?

I implemented the algorithm descussed here in the 1 D case (when $$g$$ is one variable function). Using a relaxation parameter the algorithm became very fast and converges in few iterations (n = 5, 6, , 8, ...) for smooth functions. After this I tried to implement it for 2 D functions, which is really what I need, but It became very slow. I think the problem is in the double integration, which is very slow. I tried to use the Aitken method to accelerate the process, but it doesn't work. Also Alex here used a Gauss quadrature formula to perform the integration. but it did not work well for me.

I want to know exactly why the algorithm is so slow in the 2 D case. Maybe someone used such an algorithm previously and can give me an idea on how to improve it.

The basic code is as follows (may be the one from Alex is better)

f[x_, y_] := x*(1 - x); (* The exact source term to be constructed *)

nsoleq = NDSolveValue[{D[u[t, x, y], t] ==
D[u[t, x, y], x, x] + D[u[t, x, y], y, y] + f[x, y],
u[0, x, y] == 0, u[t, x, 0] == 0, u[t, 0, y] == 0,
u[t, x, 1] == 0, u[t, 1, y] == 0},
u, {t, 0, 1}, {x, 0, 1}, {y, 0,
1}]; (* nsoleq[1,x,y] is the observation used in construction of \
the source term *)

(*The iteration is as follows g[0]=0 for example,
g[i]=g[i-1]+a[i]*p[i]
p[i] is the gradient descente of the Conjugate gradient method (we \
need to solve to pdes to calculate it)
a[i] is a relaxation parameter
*)
g[0][x_, y_] := 0; n = 20; (* Initialization of the iteration *)
Do[nsol[i] =
NDSolveValue[{D[u[t, x, y], t] ==
D[u[t, x, y], x, x] + D[u[t, x, y], y, y] + g[i - 1][x, y],
u[0, x, y] == 0, u[t, x, 0] == 0, u[t, 0, y] == 0,
u[t, x, 1] == 0, u[t, 1, y] == 0},
u, {t, 0, 1}, {x, 0, 1}, {y, 0, 1}]; (*
Solve the direct equation *)
nasol[i] =
NDSolveValue[{D[v[t, x, y], t] ==
D[v[t, x, y], x, x] + D[v[t, x, y], y, y],
v[0, x, y] == nsol[i][1, x, y] - nsoleq[1, x, y],
v[t, x, 0] == 0, v[t, 0, y] == 0, v[t, x, 1] == 0,
v[t, 1, y] == 0}, v, {t, 0, 1}, {x, 0, 1}, {y, 0, 1},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 5*15 + 1, "MaxPoints" -> 5*15 + 1,
"DifferenceOrder" -> Automatic}}];
(* Solve the adjoint equation *)
(*p[i]=N[Integrate[nasol[i][t,x,y],{t,0,1}],10]+0.000001*g[i-1][x,
y]; *)
b = N[Integrate[nasol[i][t, x, y], {t, 0, 1}], 10] +
0.000001*g[i - 1][x, y];
p[i] =
Interpolation[
Flatten[Table[{{x, y}, b}, {x, 0, 1, .1}, {y, 0, 1, .1}], 1]];
nsol1[i] =
NDSolveValue[{D[u[t, x, y], t] ==
D[u[t, x, y], x, x] + D[u[t, x, y], y, y] + p[i][x, y],
u[0, x, y] == 0, u[t, x, 0] == 0, u[t, 0, y] == 0,
u[t, x, 1] == 0, u[t, 1, y] == 0},
u, {t, 0, 1}, {x, 0, 1}, {y, 0,
1}]; (*This is for calculating the parameter of iterarion a[
i] for accelerating the convergence *)
a[i] = (NIntegrate[(p[i][x, y])^2, {x, 0, 1}, {y, 0, 1}])/(N[
NIntegrate[(nsol1[i][1, x, y])^2, {x, 0, 1}, {y, 0, 1}],
10]);(*The parameter of iterarion,
In this division we have many problems with NIntegrate,
especially when we increase the number of iteration *)
g[i] =
Interpolation[
Table[{{x, y}, g[i - 1][x, y] - a[i]*(p[i][x, y])}, {x, 0,
1, .1}, {y, 0, 1, .1}]~Flatten~1]; (*Iteration *)
, {i, 1, n}]; // AbsoluteTiming


{128.52590416078385, Null}

 {Plot3D[f[x, y], {x, 0, 1}, {y, 0, 1}, PlotLabel -> "Exact Solution"],
Plot3D[g[n][x, y], {x, 0, 1}, {y, 0, 1},
PlotLabel -> "Constructed Solution", PlotRange -> All]}


• The performance loss is due to the fact that you use interpolation function all the time. I'd suggest to get in touch with low level FEM programming and to use the coefficient vectors, and mass and stiffness matrix directly. Then the integrations in the definition of a[i] boild down to mere matrix-vector multiplications of the form p.mass.p`. – Henrik Schumacher Mar 9 at 6:55
• @Henrik, You're right. In fact I want to use this approach but I took the easy way using MOL methods (instead of FEM) and direct MMA methods since I just began MMA programming and I have a more complicated equations (not this one in the post). – S. Maths Mar 9 at 9:29
• @Henrik, van you refer to any example in this forum which use low level FEM ? Thank you. – S. Maths Mar 12 at 23:53
• The first code block of this pos shows how to obtain mass and stiffness matrix. You can find more details in the Finite Element Method User Guide. – Henrik Schumacher Mar 13 at 6:43
• You're welcome. – Henrik Schumacher Mar 13 at 13:01