So for a change, we're going to solve the one-dimensional elastic wave equation which is a differential equation in u for u is the displacement. On the left-hand side, we have the second time derivative and the density. On the right-hand side, we have mu which is the shear modulus multiplying the first derivative of u and again this whole term has to be another first derivative plus the source term. Now, what happens if we apply the Chebyshev method to such an equation? Well, on the left-hand side as you see here. Again, we'll make use of the standard finite difference operators to extrapolate the system in time. On the right-hand side to numerically solve the first derivatives, we now make matrix-vector operations and this will return the derivative and that's going to be plugged into that equation, that algorithm then to extrapolate the system in time, and again remember, we'll have an exact derivative so the only error that we commit is coming from the finite difference approximation to the second time derivative on the left-hand side. But to look at this how this works in practice, we now go to our Jupyter Notebooks. So, we apply the Chebyshev method to the one-dimensional elastic wave equation which is again an equation that has first derivative which we have to apply sequentially twice to solve the entire system. Again, remember also this time we're not calculating the derivative in the spectral domain, but we initialize a differentiation matrix and the numerical derivative is then calculated as a matrix-vector multiplication. So, let's go to the Jupyter Notebook and see the implementation of that algorithm. So, we again see here the one-dimensional elastic wave equation. The left-hand side has the second time derivative that we solve again with the standard finite difference operator a three-point operator that allows us to extrapolate the displacement at point j plus one. This is actually the future point and spatial point i. On the right-hand side, we have the divergence of the stress which is mu times the space derivative of the displacement plus a force term. Now, how do we do that? So, we have to initialize the differentiation matrix which we denote by D i j where this algorithm here, the corners actually get a special treatment, and then the formulation here is that we calculate the derivative of vector u, the displacement for the ui, as D i j u j this is an implicit summation over j that will return us the first derivative at each grid point. Now, we have a function or subroutine called get cheby matrix that initializes that matrix, so that's relatively straightforward to initialize. But let's also have a first look at this graphically. So, here is a figure of the differentiation matrix, and the important point is it's not banded or it's actually an entirely full matrix. This would be the case also in the pseudospectral Fourier implementation, the equivalent differentiation matrix in the Fourier case would actually also be full. You also see that the values actually are quite large at the edges here. This is, in fact, a consequence of the increasing grid density or the decreasing grid distance at the boundaries. We will see this later in the simulation. We skipped the initialization, you can have a closer look at this yourself and focus on the actual extrapolation scheme. So, this again is the major part of our simulation the extrapolation, and you see we are calculating this with a dot product. This is a function of NumPy applied to the pressure field that is in MPP.t means actually this is a transpose of the vector p containing the actual pressure field, and the rest is a very simple straightforward extrapolation combined with adding the source at the appropriate grid points. Actually, here we have a spatial source time function called SG that is spreading basically, a Gaussian function using a Gaussian function, the source across a couple of grid points you can- I invite you to look at this in more detail. But actually, let's just run it and see what happens. So, for a reason, this time we show the snapshots which is like you can see this looks like a string motion, not like a line, but with the actual grid points, and the reason is that that allows you to appreciate that because of the specific discretization using Chebyshev polynomials, the grid density is actually much higher at the edges than in the center, and in our case here we have a homogeneous medium. What does that mean? That actually means, at the boundaries, you sample the wave or one wavelength with many more points. So that's an overkill. You actually would not need this. So, that also means because the dx is the greatest sense is much smaller at the boundaries, you actually have to lower your timestep equivalently. So, that's kind of a disadvantage. There is one advantage and that we've not talked about this much, it actually allows boundary conditions to be implemented more accurately than compared to the Fourier methods. So that was one of the major advantages, for example, a free surface boundary condition can be implemented quite elegantly with the Chebyshev method. So finally let me say again that this limited area type calculation, it's also possible with different basis functions like Lagrange polynomials, will play a very important role later when we talk about the so-called spectral element method. But for that, we will first discuss the fundamental concepts of finite element methods in the next week.
