we have to calculate integrals like the ones that you see here. For example, integrals from minus one to one now inside your element, that contain functions of rho, density of psi multiplying the Lagrange polynomials, LI, LJ Or always functions of xi and the Jacobian. We need to calculate integrals over the stiffness values defined inside the element psi, for psi from minus one to one multiplying the derivative, the space derivatives of the Lagrange polynomials, with the terms containing the inverse Jacobian. Now it's very important that, I mentioned this before. The density and mu are allowed to vary inside the element. In this formulation, actually originally, they would be continuous functions inside the element and we will not be able to solve that numerically. In practice, what will happen, actually, we will say we know the parameters of our partial differential equation, density and mu. Only at the co-location points that we're going to use and these are the Gauss-Lobatto-Legendre points. Now a very important step, at this point, is actually we're going to use the same co-location points for our numerical integration scheme, so that's the Gauss-Lobatto-Legendre integration. And that's one of the key elements, now for this specific formulation of the spectral element method and that's the next step. So numerical integration could be the topic for an entire week. But we will focus here, basically, on just the specific method that we will use here for our scheme. Now maybe you are familiar with Gauss quadrature in which case, if we want to calculate an integral over a function, in this case, defined from minus one to one, f of psi, d psi. We approximate this by a sum over, now we're replacing the function by a polynomial approximation. And by that, actually this whole integral can be turned into a sum, which is of course very convenient and much better to do on a computer. The sum will be from i equals one to n plus one and will again be the order of the polynomial. And we have weights wi multiplying the functional value at the points psi i. Now in classic Gauss quadrature, the points actually are inside the domain and the specific aspect that the Gauss-Lobatto-Legendre polynomials or co-location policies actually that we include the boundary points, which in our case are actually the boundaries of our elements. That's in principle the overall scheme of this, which again, as I said before, by doing that, we will end up with a very compact scheme calculating the mass matrix leading to its diagonal form, that will be later. First let's make an example, and see how this numerical integration works. How well does this numerical integration work? Let's make an example. Let's take a simple function that we can analytically integrate to compare with the exact result. Now, here's an example. It's simply a sum over sum sine functions, that we say pi divided by some constant, ai, multiplying psi plus ai as a shift. Here's an example of some values that you can use. You can play around with this to find an appropriate function. Now the function is shown here. And now let's calculate the numerical integration using this GLL procedure and compare with analytical result. An example is given here for order three. The thick line is the actual function and the thin line, thin solid line is the polynomial approximation that is actually then evaluated. And you can see that by this approximation goes exactly through the co-location points. Now in this case and the error is given at the top, the error is actually very large. And that would not be an acceptable value as an approximation to the integral. Now if we increase the polynomial order, as you can see here, the value of the final integral now is already very exact and the amazing thing is the polynomial approximation is actually still not looking much like the original function. Still the integral is actually very well approximated. And you can kind of guess that this is the case by looking at the area that's the kind of the light shaded area, is the difference between the exact function and the numerical approximation, the polynomial approximation. And you can see that's kind of a zero-mean function if you integrate this difference up, you actually definitely end up with a very small value which indicates that the integral may actually be very well calculated. So that's again, a very powerful scheme as a numerical scheme to integrate a function where we do not know the analytical solution. Now you could ask the following question, but where does he get the actual co-location points from and these weights, wi? Well, again, we could make a longer story about this, but it's very easy to get the tables that you can find anywhere. And the table is shown here for up to a certain polynomial order so here the actual psi values are given and also that the polynomial weights that you have to employ to calculate numerical integrations. Now, basically, that already concludes the toolbox that we need to calculate our system matrices. We have an interpolation scheme based on Lagrange interpolation. We have a numerical integration scheme that actually works on the same co-location points. And there's one thing that's missing, and that's we have to find a way of calculating the derivatives of the Lagrange polynomials at all co-location points.
