3
$\begingroup$

I am learning finite element method(galerkin method) for solving ode/pde.

when searching this topic, I often see examples using the hat functionUnitTriangle[x] as the basis function of the galerkin approximation.

I understand that the Galekin method is a method of expressing the objective function by the sum of a basis function and a coefficient, and solving the algebraic equation that is the result of integrating the residual with each basis function.

something like, f is approximate function(solution)
u is true(target) function approximated by f
R is residual for example,$R[f(x)]=(u(x)-f(x))^2$

then

$f=a_1 \phi_1+a_2 \phi_2+a_3 \phi_3+....+a_n \phi_n$
where $\int_{R}\phi_i \phi_j =0$
if
$i\neq j$
and

for all j=1,2,...,n
$\int_{R}R[f(x)]\phi_j=0$

however,this is contrary to my intuition because linear combination of hat function that can be used for galerkin method doesn't provide practical approximation for the target function.

the below is example with 10 nodes which can't be used for galerkin method.

(*Hat function*)
kernel[j_] := UnitTriangle[x - j]
(*candicate solution*)
f = Total@Table[c[j]*kernel[j], {j, -10, 10, 1}];
(*L2 norm between target function and candicate solution*)
L2norm = Total[Power[Table[Sin[j] - (f /. x -> j), {j, -10, 10, 1}],
    2
    ]
   ];
sol = Last@NMinimize[L2norm, Table[c[j], {j, -10, 10, 1}]];
Plot[{Sin[x], f /. sol}, {x, -10, 10}]

enter image description here

in the above,inner product of basis function doesn't 0.

 Table[
  (*inner product of basis function*)
  NIntegrate[kernel[j]*kernel[j + 1], {x, -10, 10}
   ],
  {j, -10, 9}
  ]

{0.166667, 0.166667, 0.166667, 0.166667, 0.166667, 0.166667, \ 0.166667, 0.166667, 0.166667, 0.166667, 0.166667, 0.166667, 0.166667, \ 0.166667, 0.166667, 0.166667, 0.166667, 0.166667, 0.166667, 0.166667}

next is example which can be used for galerkin method where inner product of the basis functions are always 0.

(*Hat function*)
kernel[j_] := UnitTriangle[x - j]
(*candicate solution*)
f = Total@Table[c[j]*kernel[j], {j, -10, 10, 2}];
(*L2 norm between target function and candicate solution*)
L2norm = Total[Power[Table[Sin[j] - (f /. x -> j), {j, -10, 10, 1}],
    2
    ]
   ];
sol = Last@NMinimize[L2norm, Table[c[j], {j, -10, 10, 2}]];
Plot[{Sin[x], f /. sol}, {x, -10, 10}]
Table[
 (*inner product of basis function*)
 NIntegrate[kernel[j]*kernel[j + 2], {x, -10, 10}
  ],
 {j, -10, 8}
 ]

enter image description here

{0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., \ 0., 0.}

it's clear something wrong...

$\endgroup$
  • $\begingroup$ I don't get your question. Of course, the basis of shape functions need not be a orthonormal basis, neither with respect to the $L^2$-scalar product nor with respect to the $H^1$-scalar product. I doubt that anybody ever required this. This is why both the mass matrix (i.e., the Gram matrix of the $L^2$-scalar product) and the stiffness matrix (i.e., the Gram matrix of the $H^1$-scalar product) have to be computed for FEM. $\endgroup$ – Henrik Schumacher May 31 at 7:07
  • $\begingroup$ Admittedly, often a lumped mass matrix is used because it is a diagonal matrix (and thus easier to invert) and because that basically does not reduce accuracy of the method. In this case, Cea's lemma is not enough to prove convergence; one has to apply a Strang lemma in order to keep track of the consistency error due to ineaxt quadrature. $\endgroup$ – Henrik Schumacher May 31 at 7:10
  • $\begingroup$ @HenrikSchumacher So, In FEM implemented in Mathematica,basis function is not always orthonormal? $\endgroup$ – pidenet May 31 at 7:14
  • $\begingroup$ Actutally, they are never orthonormal. With very high order bases, people try to make the basis functions as orthonormal as possible (at least within the local space of an element) in order to obtain sparser matrices. However, Mathematica uses basis functions of at most order 2 and orthogonalization techniques just do not pay off for such low degree - so I think. @user21 may be able to give some more detail; he is the developer of the Mathematica's FEM framework. $\endgroup$ – Henrik Schumacher May 31 at 7:14
  • $\begingroup$ your comment immediately solved my biggest headache this year...thanks! $\endgroup$ – pidenet May 31 at 7:16
2
$\begingroup$

Turning my comments into an answer.

Of course, the basis of shape functions need not be a orthonormal basis, neither with respect to the $L^2$-scalar product nor with respect to the $H^1$-scalar product. This is why both the mass matrix (i.e., the Gram matrix of the $L^2$-scalar product) and the stiffness matrix (i.e., the Gram matrix of the $H^1$-scalar product) have to be computed for FEM. Actually, one needs only the mass and the stiffness matrix; they a assembled from elementwise contributions (the so-called local matrices) and the nodal functions themselves are never evaluated explicitely (or only at a few quadrature points).

The set of nodal basis functions is never orthonormal. With very high order bases, people try to make the basis functions as orthonormal as possible (at least within the local space of an element) in order to obtain sparser matrices. This is however quite an algebraic fuzz. Mathematica uses basis functions of order at most 2, and orthogonalization techniques just do not pay off for such low degree - so I think. @user21 may be able to give some more detail; he is the developer of the Mathematica's FEM framework.

What's important is that the nodal shape functions of two nodes that are far apart are orthogonal to each other; otherwise, mass and stiffness matrix would be dense matrices and thus too large to be stored in memory.

With the finite elements in Mathematica (which are $C^0$), two nodes are "far apart" if they do not belong to the same finite element. For $C^k$ elements with $k>0$, stencils of a node are larger (e.g. the 2-ring of a node for $C^1$). Thus, the higher the $k$ and the higher the order of the bases function, the less sparse the mass and stiffness matrix are. This is why one usually tries to use as low degree elements as possible, in particular when the solutions of the PDE are too rough to enjoy faster convergence rates with high order basis functions.

Finally, I would like to mention that a some contexts, a lumped mass matrix is used because it is a diagonal matrix (and thus easier to invert) and because that basically does not reduce accuracy of the method (as long as the solution of the PDE is sufficiently smooth). In this case, Cea's lemma is not enough to prove convergence; one has to apply a Strang lemma in order to keep track of the consistency error due to ineaxt quadrature.

$\endgroup$

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.