[MUSIC] So before proceeding to the final algorithm, let's again, look where where we come from. Here's again the weak form in discrete version of the elastic wave equation where we have sums over integrals containing the basis function themselves and also sums over an integral containing the derivative of the basis functions. Now we introduced the mass matrix and the stiffness matrix. We have the vector u of the unknown displacements. And we also have a vector describing the source term. Now, we can get rid of the sum signs of course by using implicit matrix vector notation, and basically arrive at the following very simple equation now, which is the second derivative of the unknown displacement in vector form, multiplying the mass matrix. Plus we have the displacement vector multiplying the stiffness matrix, and on the right-hand side we have the force terms. Now, the only slight modification we do is actually we transpose the mass and the stiffness matrices to arrive to the final form as it's given here. So basically, this is the solution, the equation that we're now trying to solve, but there is one problem, we have the second derivative of the displacement values, the vector of displacement, what are we going to do with this? So we have to first solve that problem, it's going to be the next step. So our final algorithm, you can see here, contains a term where we multiply the transpose of the mass matrix with the second derivatives of the displacement vector. Now, we simply apply the finite difference scheme to the second time derivative of the displacement vector, and you can see that here. So we take basically u(t + dt)- 2u(t) + u(t- dt) divided by dt squared. And plugging that into the equation, we can basically now as we've done several times, with basically all our schemes that we've encountered so far, we can take out the value of the displacement field in the future at (t + dt). And express everything in terms of displacement values that we know at u(t), which is the present time, or u(t- dt) which is basically the past, one time step in the past. So the final algorithm that we're actually going to solve, and this is more or less than exactly written in that form in our python code, you can see here we have u(t + dt), given as a function. Now, as you can see here this is the inverse of the mass matrix, plus all the other term that you see on the right hand side. It looks very, very similar to the finite difference type algorithm and actually we will show before we go to the Python code that formally you can express a classic final difference algorithm in exactly the same way, introducing a differentiation matrix for the finite difference scheme. Before proceeding to implement the algorithm, we have to actually find a way of solving and calculating the mass matrix and the stiffness matrices. For that, as mentioned before, we make use of an essential ingredient of basically all finite element type algorithms, we move to a local coordinate definition. And that's given here, again, like in the static case, we replace x by a coordinate xi, which is defined as x- xi, and we now become more flexible actually. We'll allow the size of the elements to vary and define the size of element i as hi, which is then of course the xi + 1, the right boundary minus xi, the left boundary. And with this definition, we can now write down the basis functions as given here now using the xi coordinate. And it also allow us to actually calculate the derivatives of the basis functions remember, we need those for the stiffness matrix. We can write it down as you see here. And an example which is more powerful in illustrating this concept, is shown here in the figure where we show both the actual basis functions, remember they are linear basis functions. So inside one element, there are just basically a straight line, connecting one boundary with the other boundary. And also we have the derivatives shown in the same graph. The derivatives of course then are constant inside each element. So with this definition of the basis functions, analytically now we know the derivatives, analytically we can now proceed in calculating explicitly the system matrices that we need for the algorithm. An important point, we can now basically calculate mass and stiffness matrices, and we do that basically as an initialization step. So we do this before starting the time extrapolation that we've just discussed in connection with the finite difference approximation to the vector of displacements. And now let's proceed and explicitly calculate some of the elements of mass and stiffness matrices. So let's start with the mass matrix. Remember the mass matrix are integrals over the density, multiplying the basis functions phi i, phi-j. Now, if we would allow an arbitrary function of space for the density inside the integral, you would actually have to apply a numerical integration to solve these integrals. So we make life substantially more simple by allowing the density to be constant inside the element. That's actually not such a strong restriction. We will relax this later when we use other methods like the spectral element method. So we take out the density of the integral, and now we're left with integrals over the basis functions phi i, phi j only. So the next one is we move from a global coordinate system to a local coordinate system, and if you recall the picture showing the basis functions and the derivatives. If we make an example to calculate an element of the mass matrix in the diagonal, Mii, it's easy to see that we actually have to break up the integral to one integral for the basically i -1 element, and plus an integral for the ith element, and you can see this here. Now, to solve these integrals is really very, very easy and you end up with a term that contains a product of density and element size for these various elements. So that's a diagonal element of the mass matrix. We live it to an exercise to actually calculate the non-diagonal elements also quite simple, but it's a good exercise to do. Now, this basically concludes the mass matrix but let's have a quick look at what the overall mass matrix looks like assuming for the moment constant element size. And you can see this here up to a factor, we have a banded mass matrix, so a very simple structure. And remember we have to calculate the inverse of this matrix before starting the time extrapolation algorithm using to find that different scheme that we discussed in the last lecture. What about the stiffness matrix? Well, it's actually very similar to the static case, except that at this point we allow the size of the element to vary. So we have integrals over mu times the space derivatives of the basis functions. So we move again also to the local coordinate system xi. And to make life easy, just like in the mass matrix case, we actually allow the string modulus to be constant inside the elements, so we can take it out of the integral. And that allows us to analytically solve these integrals, which is very, very handy. So we show an example here for a diagonal element of the stiffness matrix and as you can see here, we again have to split it up into two integrals, one of the elements i- 1 and 1 over element i. And we end up with terms containing the shear modulus values of these elements divided in this case by the element size. So again, I'll leave that as an exercise to evaluate the off diagonal elements, very simple but this is a really nice way of getting into how these matrices are calculated. And then we basically have all these ingredients now to solve our system, and basically solve the one dimensional elastic wave equation using a python algorithm. Before we move to the final algorithm, something else, we would like to compare with the finite-difference solution. It's always nice if we can compare a new scheme to something we know very well which is the finite difference solution. It turns out and we have seen this in when we discussed the Chebychef method, Chebychef pseudospectral method, a derivative can also be written as a matrix vector operation. So basically here on the left-hand side, you would have the vector of displacement, the second space derivatives, and that's equal to a multiplication of a differentiation matrix D multiplying the vector of displacement u. Now, what does this differentiation matrix has to look like in order to do this operation?
