So, back to physics before we start solving these equations numerically. So, the equation that we're going to use to illustrate the various numerical methods is the wave equation. Now, just a few words on what to expect in terms of solutions. So, if I write the equation again here which is the second time derivative of p, function of x and t, is equal to c square, phase velocity square. Multiplying the second space derivative of the pressure function of x of t plus a source term also depending on x and t. Now, there are two situations where we can actually find an analytical solution. So, if c is not a function of x, we can and actually now we also avoid the source term, so no source term. We have with an initial condition P at t equals zero equals P_0. Also the gradient of P being zero, we can find analytical solutions, we will derive that later and show you the result. The second situation is, if the source term s of x, t actually is a delta function in x and t. So, as you probably know that means we will get the impulse response of our partial differential equation in terms of a linear system, we will talk about linear systems again at the end of this week, and also as a very important concept in connection to solving partial differential equations. So, the impulse response is basically the Green's function. So, if we have on the right-hand side a source term in terms of delta functions in x and t, then we can basically replace the P with a capital G which is the Green's function, the solution to an impulse response. We're not deriving this here but we will make use of these analytical results many times. Why? Because in this situation where we have a homogeneous medium so, the velocity does not depend on space, we actually have analytical solutions to compare our numerical solutions with and that's very powerful. It's extremely useful to test whether what you do mathematically, if you develop an algorithm and if you implement it as a computer code whether you actually have a correct result. We'll also make use of a slightly different form of the wave equation, this is for elastic waves. You see here some strings vibrating, it's basically a description of waves on a guitar string for example. That's described with the letter u, which is a displacement. In this case, in the y direction, if you propagate in x. Let say this is x, then you have a particle motion perpendicular to the propagation direction, that's a transversely polarized wave. So, the partial differential equation describing this is the second time derivative of u_y multiplied by Rho, is equal to the first space derivative with respect to x multiplying the open brackets Mu, which is the shear modulus which is a function of x. Also, the density is a function of x multiplying the first space derivative of u plus some source term that acts in the y direction. Now, something interesting happens if you again assume a homogeneous medium. So, density and the shear modulus Mu does not depend on space, then you can take the Mu out of the bracket. If you take Rho to the other side of the equation you get Mu divided by Rho. You might notice, the shear velocity actually is, let's call this c, is equal to the square root of Mu divided by Rho. So, you end up with an equation second time derivative with respect to time of Mu is equal to c square times second time derivative of u plus a source term and you see, mathematically that's actually identical to the wave equation that we had for acoustic waves. So, on the other hand, when the medium is heterogeneous, this equation is a little more complicated than the acoustic wave equation as you've seen because you also have a derivative of something multiplying the shear modulus. We will get into that later. That actually has an impact on the specific solution algorithm. To understand the actual original real partial differential equation, but also to understand the properties of the numerical algorithms that we developed, we often make use of so-called plane wave descriptions. You've encountered this before, but just one more word on this. So, a plane wave is characterized by a wave number k, which is basically a vector that points in the direction of propagation as you can see here. So, k is a vector kx, ky, kz components. So, if you think of a pressure, a plane pressure wave, this would be characterized as p is equal to p_0 which is an amplitude multiplying for example, sine kx, where k is a vector, x is a position vector minus Omega t. It could also be a cosine replacing the sign or a combination of both. We also often make use of the exponential form , in which case we would write p equals p_0, e to the ikx minus Omega t after i open bracket of course. Again, k is the vector, the wave number vector, x is the position vector. The same applies to the elastic form. In this case for example, the description of an elastic transverse wave in one dimension Let's make an example, like this would be u_y equal to A_0. There will be an amplitude in the Y direction, e to the i k_x, x minus Omega t and kx, then would be the wave number in the x-direction. So, all the time here we have the relation C equals Omega divided by k. If it's a vector, then it's actually the modulus of k. If we have a one-dimensional wave propagation, C equals Omega divided by k_x if we propagate in x-direction. Our linear wave equation actually has a very interesting property, it's symmetric in time. Now, I would like to illustrate that with an example what that means. You see here one-dimensional example of a complicated model, so, it's the structural heterogeneity here that the velocity model as a function of space is very complicated, it's actually random. So, let's assume that we have two points, A and B. We start with injecting a source at A recording at B, then we inject a source at B recording at A. Would you immediately think that this gives exactly the same result? I don't think it's so intuitive, but actually that's the result the consequence of the fact that the wave equation is symmetric in time. So, this is a numerical example, and so you can redo this later once we start developing our own algorithms. So, it turns out you actually get the same result. Another aspect of this is the so-called Time reversal. Also, behind me you will see a wave field that's propagating out from a source point, and it's actually recorded at your circular array of receivers. Now, you can actually take these recordings at each of those receivers in the ring, you'd flip them around in time and re-inject them as sources at those points and you can see what happens. Actually everything is propagating back to the source and focuses on the individual source point. To some extent, you can call this reverse acoustics and actually that's something that's for example used in medicine to focus acoustic energy in sonic applications. But, the mathematical concept of this time-reversal reciprocity is actually also extremely important in developing tomography algorithms, meaning recovering the internal structure where you have recorded acoustic waves, for example in oil exploration problems So, again, these are fundamental concepts of the wave equation it's real fun playing around with these things using numerical algorithms. Let's talk a little more about structural heterogeneities. What I mean by that is basically the variation of the parameters of the partial differential equation. In our case, this can be the acoustic velocity. Basically, the speed of the waves as a function of space or it could be in the elastic case, the distribution of density and the shear modulus. Now, here's a graph that shows some examples of different types of heterogeneity. If you think of the Earth, this could be a layered model, it could be for example, a model with internal curved interfaces, where the material changes abruptly, or you might be interested in understanding the scattering of a single point, where the material parameters are changing. There are also different degrees of smoothness of heterogeneity. So, for example, there might be a relatively smooth model of heterogeneity's, and in that case, there might be other more efficient techniques than directly solving the wave equation using numerical methods. So, what I'm thinking about is for example, ray theory or finite frequency ray theory. In the most general case and that is given here on the bottom right, this would be basically a medium with arbitrary heterogeneities, random heterogeneities for example, and there's just no way around using numerical methods. Again, I want to recall that's probably the reason why you're here because you want to solve a problem where there actually is a certain degree of structural heterogeneity of the partial differential equation we're going to solve, the parameters are actually varying rapidly in space. There's another concept I would like to explain to you in connection with solutions to the exact equation and our numerical approximation. So, remember, we would like to find solutions to the wave equation. Let us take the acoustic wave equation and we look for pressure as a function of time. When a homogeneous medium, and we know for these kind of situations, we have actually a unit function. We have an analytical solution. Now, without showing equations here, I just want to explain the concepts. Actually, the analytical solution to the problem of delta function, is a Heaviside function, which is the integral of a delta function. So, actually what you expect, basically for an impulse at some distance is just simply a step function that's given here. Now, would we ever be able to simulate something like this correctly with a numerical solver? The answer is no. Why? Because the Heaviside function or the delta function as you probably know contains infinite frequencies. We will never be able to simulate on a computer in the discrete world infinite frequencies. So, we'll always have a band limited solution to the Green's function. Now, I would like to show you something interesting that happens. But, if you have a numerical solver and you inject such a delta function. Actually, you can do this. We will see this later, and I just want to show you the results here. So, the result is something like you see here, and it looks weird. It has a step, by then it has a lot of oscillations. That's called numerical dispersion, we'll understand this later. But, does that mean we have to throw away this solution, we cannot use this? The answer is no. Actually, as you know, there is the so-called, convolution theorem that tells us, that if you want to obtain the solution of a system for an arbitrary source time function, you actually do a convolution of the Green's function with the source time function. You see, this funny symbol in between G and the source time function, which denotes the convolution. There is also, remember the convolution theorem that tells us that if there's something that is defined a convolution in the time domain, then in the frequency domain, it's a multiplication. Now, what happens if you obtain such a numerical solution that obviously seems to be wrong? Now, let's assume you would like to know the response of the system to let's say, the first derivative of a Gaussian. This is our source time function. The analytical solution is simply, the convolution of the Green's function with this source time function and that's what you obtain, is actually the integral of the first derivative of a Gaussian. So, it will be a Gaussian function. Now, what happens if you do exactly the same operation and not with the exact Green's function, but with the numerical Green's function, which often you see seems to be wrong. Surprise? We actually obtain the correct solution because you're filtering out the part of the solution that is incorrect, you still obtain in this case, for these parameters, the correct solution of this problem. Now, I hope I haven't confused you, that might seem academic, but it actually has a lot of potential applications. Because that means, if for example, you want to understand wave propagation in a very strongly heterogeneous, maybe a random medium, you can actually inject a delta function in your numerical algorithm, you get basically completely wrong numerical solution, but you can filter it later. If you make sure that your algorithm is correctly modelling the frequency domain of your source wavelet, you will actually obtain the correct result. So, that's a consequence of the fact that not only the partial differential equation itself, but also our numerical solution can be treated as a linear system with all the consequences concerning Green's functions and convolution with potential source time functions. So, a very powerful concept.