So, we now have basically everything to turn this into a computer program, but what to compare with. I want to show you an interesting, other approach using the finite-difference method to solve such a static problem. So, let's first look again at the original equation. So, we have minus Mu, multiplying the Laplace operator, equaling a force term on the right-hand side. Now, what happens if we simply replace that Laplace operator, so the second derivative of u with a finite difference approximation? Of course, we use a three-point, simple three-point operator here. So, that will be ui plus one, minus 2ui, plus ui minus one, divided by h squared. Let's use h here as the distance between two elements and we keep it constant. So actually, we're now solving this numerical approximation for ui. What happens as you see here, you end up on the right-hand side with an average using the adjacent values, divided by two, that's an average basically, plus a forcing term. Now, the way one can solve for u is now, you start with an initial guess, for u, it can be zero, and then you basically iterate it. So, you start using the update for u each time, just a little bit like a time evolution, but in this case, it's actually an iterative approach. Here's a very simple piece of Python code where you can see how this is implemented. Looks very simple, just a couple of lines where basically again, you're iterating, and of course, you have to make sure that you iterate long enough, that you converge to the right solution. How do you know that? Well, you can compare with the finite element method, and we now know how to do this. These two approaches, we now compare inside a Python code in our Jupyter Notebooks and look how they perform. So, we go back to our Jupyter Notebooks, and I hope you have also as much fun as I have developing and using these Jupyter Notebooks. So, we'll solve our first finite element type problem. First, the static problem. For the first time, we will encounter the fact that we need actually matrix inversion to solve our system. Let's have a look at our code. We start again with the original equation minus Mu times the second space derivative of u, is equal to the forcing term. We were able, using the finite element type algorithm, the Galerkin principle to write this as a matrix, as a linear system of equation using matrix vector notation. So, we're solving for u, we need to calculate the inverse of our stiffness matrix K or better the transpose of K, and applying it to the forcing f, forcing vector f. The stiffness matrix again, as you can see here has a banded structure, the numbers here corresponding to the finite difference, the second finite difference approximation, apart from a factor of minus one. It's also not divided by h squared, as it is in the finite difference case by h. Now, we will compare this with the relaxation method that we just described in the course, and we update it in an iterative way to obtain the solution and that allows us to compare with the finite element result. Now, we initialize the spatial domain with 20 grid points. Remember, these grid points are actually the element boundaries. So, we will in fact, we'll have nx minus one elements. So, we would have 19 elements. We initialize with the same shape, the solution vector u, and the source vector f. In this very simple example, we keep the stiffness mu, the shear modulus mu constant at one, for the entire domain, which does not depend on space here. We initialize space, which is defined here between zero and one with nx points, and that also allows us to calculate the element size simply by taking the difference between two element boundaries. The stiffness matrix and that's a classic thing that we have to do in preparation for any finite element solution, is initiated here using for loops. You can see that depending on where we are inside our matrix, the diagonal or non-diagonal elements, we can simply initialize it with the values we calculated previously in the lecture. We put a source term in space at three quarter of the entire domain. So, three times nx over four. Also just to show how this works with the boundary conditions, we impose a boundary condition on the left side upon zero to 0.15 and on the right side, we'll have another boundary value u, 0.05. Actually, as you can see here, we simply turn this into an additional term for the force vector at the boundaries. Now, the finite element solution is now basically written in one line, and that line, as you can see here, which is the Python implementation of u equals inverse matrix K transpose, times forcing vector, and we are basically done, and we can directly plot that. But before that, we compare it with the finite difference solution, the relaxation method as we've described before. You've already seen part of the Python code in the course, and it's again shown here, we just initialize u as zero, and then update it iteratively until we converge. Now, let's see what this looks like. So, we run it. Then you can see that the dashed dot line this is the finite difference approximation, and the red line is actually the finite element approximation. Very slowly, sorry, converges to the finite element solution. So, that's basically, all about the static case. We solved the problem after initialization, the stiffness matrix basically with one line, which is very elegant. The relaxation method takes longer to calculate, but again, it doesn't mean one is better or worse, it's simply a way of showing the working of two entirely different methods to come up with a solution here for the static elasticity case.
