0
$\begingroup$

I need to explain the background of the problem I'm working on, and the questions themselves are a couple of paragraphs down. I am reading this book and implementing the algorithms along (so far in chapter one, assemble the mass matrix and the load vector in 1D). As I go on, I'm writing unit tests. The one for the assembly of the load vector is failing, either because of a programming bug, or because I'm failing to understand the math. Please let me know if it's the latter.

At his point the book isn't laying out any differential equation to solve. It is just defining the mass matrix and load vector, on the mesh/interval $I$, in the context of approximation and projection. In this context the mass matrix is defined by the book as $$ M_{ij} = \int_I \phi_j \phi_i \,dx $$ and the load vector for a function $f(x)$ as $$ b_i = \int_I f \phi_i \,dx $$ where $\phi_i$ are the "hat functions," basis of the subspace of piecewise linear functions on the mesh of $I$.

So the whole point of these definitions in this chapter of the book (again, not yet laying out any differential equation to solve) was supposedly that the equation $$ \mathbf{M \xi = b} $$ holds for the coordinates $\xi_i$ of the $L^2$-projection of $f(x)$ onto the mesh.

QUESTION 1: is this so? Have I understood the book correctly so far?


On this assumption I based my unit tests. First, I understand that $$ \sum_i \xi_i \phi_i(x) \approx f(x) $$ Moreover, I understand that if I choose a mesh $I$ and a function $f$ which happens to be already piecewise linear on $I$, the approximation should become identity.

Therefore I made up the following example: $$ \mathbf{x} = \{ -10, -4, 8 \} \\ f(x) = \left| 2 x + 8 \right| $$ which gives me the following mass matrix and load vector: $$ \mathbf{M} = \begin{bmatrix} 2 & 1 & 0 \\ 1 & 6 & 2 \\ 0 & 2 & 4 \end{bmatrix} \\ \mathbf{b} = \{ 36, 0, 144 \} $$ Since the projection should equal the original function, because the latter was already in the subspace (*): $$ \mathbf{M \cdot \xi \stackrel{*}{=} M} \cdot \{ f\left(\mathbf{x}(1)\right), f\left(\mathbf{x}(2)\right), f\left(\mathbf{x}(3)\right) \} $$ However this gives $$ \mathbf{M} \cdot \{ 12, 0, 24 \} = \{ 24, 60, 96 \} \neq \mathbf{b} $$ Not even approximate.

QUESTION 2: what have I done or understood wrong?


The calculation of the load vector is done numerically. The book uses the Trapezoidal rule, and proposes the following Matlab code:

function b = LoadAssembler1D(x,f)
n = length(x)-1;
b = zeros(n+1,1);
for i = 1:n
    h = x(i+1) - x(i);
    b(i) = b(i) + f(x(i))*h/2;
    b(i+1) = b(i+1) + f(x(i+1))*h/2;
end

QUESTION 3: is this implementation correct, assuming the Trapezoidal approximate integration?


This code computes for this mesh and $f(x)$ the following load vector:

>> LoadAssembler1D( [ -10, -4, 8 ], inline('abs(2 * x + 8)', 'x') )
ans =
    36
     0
   144

where I got the discrepancy above. What's wrong? Is it the implementation, or something else?

Is my mesh perhaps so coarse that the numerical integration may be to blame for this whole error, even if the code is correct based on it? I had the idea that the Trapezoidal rule would be exact for this simple case, but now I realize that my load vector is made up of integrals of piecewise quadratic, not linear, functions...

Or is there something else wrong that I've missed?

$\endgroup$
  • $\begingroup$ The mesh is too coarse. I think you are just getting a huge integration error while integrating the second degree polynomials ($f \phi_i$) with a first degree quadrature. Normally, most FE solvers (especially for 2D, 3D problems) will assemble the matrices pulling back to a reference element, where Gauss quadrature with higher precision is used. $\endgroup$ – PierreCarre Feb 6 at 15:30
  • $\begingroup$ For instance you are considering $\int_{-4}^8 \frac{1}{12}(2x+8)(x+4) \,dx \approx \frac{12}{2}(0+24) = 144$, which is not even close to the real value. $\endgroup$ – PierreCarre Feb 6 at 15:41
  • $\begingroup$ Indeed; updated the answer... thank you! $\endgroup$ – J.P. Feb 6 at 19:38
0
$\begingroup$

You have understood correctly. The $L^2$ projection, $f_h$, is the solution of the problem $\langle f_h - f, v_h \rangle =0, \forall v_h \in V_h$. Writing $f_h$ in terms of the basis of $V_h$, $f_h = \sum_i \xi_i \phi_i$ and choosing $v_h= \phi_1, \cdots \phi_n$ you get the linear system you mentioned. So yes, the solution of $M \xi = b$ should give you the $L^2$ projection.

I think your load vector is not properly computed.

$$ b_i = \int_{-10}^8 f(x) \phi_i(x)\,dx $$

yields $b = (24, 60, 96)$ (see the note below), and the solution of the system $M \xi = b$ is $\xi = (12,0,24)$, as expected.

$\textbf{Note}$: for instance $$ b_3 = \int_{-4}^8 \frac{1}{12}(2x+8)(x+4) dx = 96. $$

$\endgroup$
  • $\begingroup$ Thank you! That's good news, but then my problem is contained somewhere else... I have edited the question to complete my story. My load vector calculation is numerical, as I said. Is the code above wrong? Where is the difference coming from, then? $\endgroup$ – J.P. Feb 6 at 15:00
  • $\begingroup$ Indeed! The only error is from the numerical integration over the super coarse mesh, as you say in the question's comment... I knew about the mesh, but the problem was that I was flying with the assumption that, since I had chosen a linear f(x), the Trapezoidal approximation would become exact. I didn't realize that what I was integrating was not just f but multiplied by the hat functions, which add another degree up to quadratic... But it's exact if I decrease one degree of f from linear to constant! Final code: github.com/XavierAP/femd/blob/master/src/femd/assembly.d#L128 $\endgroup$ – J.P. Feb 6 at 19:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.