[MUSIC] So we start again, as usual, with our one dimensional wave equation, elastic wave equation. Density multiplying the acceleration on the left hand side. We have the derivative of the stress on the right hand side plus a source term. Now, we start with the classic derivation steps for Galerkin type methods. We first integrate over the entire domain D, and also at the same time multiply the equation with an unknown test function that we call phi. We're still entirely in the continuous world. Now, we'll take the first term on the right hand side, we take on the left side. So we have basically the integral containing the acceleration minus the integral containing the shear modulus and the remaining terms on the left hand side. And to this now second term, we apply partial integration. Now, as we've seen in the finite element derivation, it turns out that we have an anti-derivative where we evaluate basically the stress at both boundaries of the domain. And if we assume it's a stress free boundary which in seismology, at least in our case, we actually usually assume, then this anti-derivative vanishes. And we're left with a positive integral, now, the minus signs turns into a positive signs. So we're left with an integral that contains the derivative of the bases function phi which we have to deal with later. So that's the first step. First thing we do now is we approximate our unknown displacement field, which in the original equation is continuous. And to get into the discrete world we replace it by an approximation u bar, which is a sum over some basic functions that we have not specified yet. So it's u i which are the coefficients, the actual values of the displacement field at some grid points inside the element. Multiplying the basis functions phi i. Now, here the sum actually goes from one to what we call capital N subscript p. Now, this capital N subscript p actually stands for the number of points for polynomial order capital N and p stands for polynomial. Because now, later on we are going to actually use polynomials, like Lagrange polynomial that go up to arbitrary order. And the number which we have the number of terms we have to sum will depend on the order of course, but it's not exactly the order, it would actually be order plus 1. So if we for example want to approximate up to order 2 we will actually have to sum 3 terms on the right inside. So with this approximation we go to the weak form of the wave equation that we have just derived. And again we take the constant value u, i out of the interval and we also take the sum sign out of the integral. And you see the now quite already lengthy equation here. But what we recognize again is the well known matrices and matrix vector form of this equation. This equation in fact they're now, again, Np equations in J and we will later see this in more detail. But what we will do now is basically identify again the matrices and vectors in that system of equation, and write it in a much denser way So with this matrix notation we can transform this animal of an equation into something that looks much, much simpler. So here we have the mass matrix M multiplying the time dependent coefficients of the displacement, u, it's basically the second derivative. Plus the stiffness matrix, capital K, multiplying the displacement vector and on the right hand inside we have the forcing, the source terms. Now, it's actually cool to look at this equation in a graphical way. And that's done here, and that's not actually, that's not an artistic view, that's an actual output of a simulation. So you have on the left hand side the mass matrix which it turns out is a diagonal matrix multiplying the vector of displacement field here or the acceleration here. Plus the stiffness matrix which is a banded matrix similar to what we had in the linear finite element method multiplying the displacement field. And on the right inside, you have a forcing term, where at some point you are injecting a force somewhere in your spatial domain. So what remains to be done here is actually our unknown field is the vector of displacement at each time step. But it's given here, still with a second time derivative and then we have to get rid of this again using the finite-difference scheme. Now, as we have done before, we now just simply have replaced this second derivative by a finite-difference operator, 3 point finite-difference stencil and by defining u new equals the vector of displacement at t plus dt. U will be the [COUGH] vector of displacement at ut, and u old will be the past, will be the vector of displacement at t minus dt. If you plug that into that equation and replace the second time derivative, we can come up with our final extrapolation scheme that basically provides us with an algorithm for the solution of this system extrapolating it in time from 0 initial conditions. And that's going to be the next step, what's remaining is of course, we have to find a way of calculating this, the mass matrix and the stiffness matrix scheme for our specific spectral-element scheme Before proceeding, let's have a quick look again at the matrices. So we have the mass matrix containing in each element integrals over terms containing the density, multiplication the basis functions, j basis functions i. So the stiffness matrix will be integrals over mu, which still is here, a continuous function, multiplying the derivatives, the space derivatives of the basis functions. Now, what is very important is for the spectral-element method, we will not take mu, hence the density, out of the integral. So what that means is basically because there is no way that for an arbitrary medium we will have an analytical solution for the integrals over lambda, and, sorry over rho and mu. It would be easy as well as the case in the linear finite element methods, to analytically solve the integral simply over the basis functions, but not if we have the parameters inside. That’s why later we’ll actually going to have to resort to numerical integration. At this point, basically, we have integrals over the entire physical domain, that’s not what we want. We need to go down to the element level to basically make all this calculations of the matrices more simple
