Implicitly restating a classical, double integral, using FreeFem.

One of the programs which I’ve been experimenting with, is called “FreeFem”. Actually, it exists both as a library, as well as a set of executables, the latter of which are meant to facilitate its use, in that scripts can be written in a loose but C++-like syntax, that will already produce interesting plots. But, plots can also be meaningless, unless the user understands how they were generated, as well as the underlying Math.

In an earlier posting, I already wrote that a basic weakness of the ~standard~ 2D integrals which FreeFem will solve for, was, that those do not correspond closely to “Double Integrals”, as I was taught those in Calculus 2.




Just to recap basic Calculus 2: If an integral is supposed to be plotted over two dimensions, then the standard way to do so is first, to declare a variable that defines one of the dimensions, and to declare a second variable that defines the second. Next, a 1-dimensional integral is defined over the inner variable, which results in a continuous function. Next, this continuous function is integrated a second time, over the outer variable, thereby resulting in a double integral.

Even though this could be extended to more dimensions than 2 or 3, there is usually little practical value in actually doing so, at least, in my limited life.

Further, the fact is clear to me that integrals exist, which go beyond this basic definition. For example, “Curl Integrals” also exist… However, for the moment I’m focusing on how to overcome some of the limitations which are imposed by FreeFem. And while working on this problem, I have found a way to force FreeFem to compute the sort of double integral which I was taught. The script below shows how I did this…


// For label definition.
int Bottom=1;
int Right=2;
int Top=3;
int Left=4;

// The triangulated domain Th is on the left side of its boundary.
mesh Th = square(10, 10, [(x*2)-1., (y*2)-1.]);
plot(Th, ps="ThRectX.eps", wait = true);

// Define a function f.
func f = x;

// The finite element space defined over Th is called Vh here.
fespace Vh(Th, P2);
Vh u, v;	// Declare u and v as piecewise P2-continuous functions.

// Get the clock in seconds.
real cpu=clock();

// Define a simple PDE...
solve SimpleXInt(u, v, solver=LU)
	= int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))	// Show me a kludge.
	- int2d(Th)(
	+ on(Bottom, u=0)   // The Dirichlet boundary condition
	+ on(Right, u=0)
	+ on(Top, u=0)
	+ on(Left, u=0);

// Plot the result...
plot(u, ps="RectX.eps", value=true, wait = true);

// Now, let's try to plot a double integral...
Vh A, Fh;

func real f2(real xx, real yy) {
//	return xx * yy;
	return xx;

func real Fx() {
	real s = 0.;
	for (int i = -5; (i * 0.2) < x; ++i) {
		s += int1d(Th, Left) ( f2(i * 0.2, y) );
	return s;

// Compute the x-integrals for reference purposes...
Fh = Fx();

solve Double(A, v)
	= int2d(Th) (
		  dy(A) * dy(v)
	- int2d(Th) (
		Fx() * v
	+ on(Bottom, A=0)   // The Dirichlet boundary condition
	+ on(Right, A=0)
	+ on(Top, A=0)
	+ on(Left, A=0);

plot(A, ps="DoubleIntX.eps", value=true, wait=true);
plot(Fh, value=true);

// Display the total computational time.
cout << "CPU time = " << (clock()-cpu) << endl;


This script actually outputs 4 plots and saves the first 3. But, for the sake of this posting, I’m going to assume that there is no reason to show the reader a plot, of a plain, rectangular mesh… The following plot shows the result, when ‘the Poisson Equation’ is simply given for a function of (x), such that:

f(x) = x




A simple linear function was given. It can be integrated once or twice, and the plot above shows how FreeFem usually does so, resulting in u(x,y).

One interesting fact about FreeFem is, that this program can be made to compute a 1-dimensional integral explicitly, i.e., without imposing any need, to solve an equation. The resulting sum (a real number) can then simply be returned with a variable, for every maximum value the function was summed to. So, what would be tempting to do next is, to integrate this variable a second time, which might result in a 2-dimensional array. However, there is a basic limitation in how FreeFem works, that would next stand in our way… (:1)

And so, in order for the second (outer) integral actually to be with respect to (y), I used the little trick in the script above. I defined the gradient of Finite Element function (A), with respect to (y) only, and stipulated that it must equal the 1-dimensional integral which was computed before, for all combinations of (x,y). Hence, this second step of my solution is similar to a degenerate Poissson equation, in that the addition of the (x) gradient has simply been dropped.

When this approach gets applied, the following plot results:




Now, there is still something which FreeFem does, which is supposed to be a feature, and which I cannot switch off. That is, to impose that at the boundaries, a certain exact value needs to exist for (A), just as was needed for (u). Documentation states, that by merely defining an equality, the (2D) derivative of a function is supposed to have with some other function, one has not fully defined that function. A boundary value must also be given, and then the resulting “Partial Differential Equation” (‘PDE’) actually defines, either the FE function (u), or the FE function (A). (:2)  Further, the way ‘int1d()‘ is implemented, also applies such boundaries.

For that reason, I can also not just shut off, anything which I feel FreeFem might be doing wrong, every time it tries to solve a PDE. What I can do is, request that this boundary occupy the left-hand side of the rectangular plot, and additionally, make sure that I begin my summation with the same value, thus resulting in a definite and not an indefinite integral every time. Also, I had suspected somewhere, that in order to resolve this issue with 2D plots, it really only needs to be resolved along one inconspicuous axis. In this case, I resolved it explicitly for the X-axis, and the behaviour of the Y-axis fell in line.

Yet, what I can do is observe, that the two plots shown above don’t match. And as long as they don’t, what FreeFem computed in its first plot, was also not a double integral.

(Updated 7/25/2021, 17h30… )

Continue reading Implicitly restating a classical, double integral, using FreeFem.