So… what is a simulation? Pt II
In the last article we were exploring the idea of how we tame the laws of physics through mathematical tricks to approximate solutions to very difficult equations. We were exploring the idea of not solving the differential equation per se, but solving a set of integral equations whose solutions are equivalent to the solution of the differential equation problem.
To recap, here’s what we had:
We were left hanging on an idea that we should approximate our solution u through a series of piecewise polynomials:
But we still don’t know what are those polynomials or what is v and how everything is connected. So let’s start from there.
A very common idea in statistics is one of regression. We take a set of discrete points and say that there is a mathematical formula that connects them. That’s my definition of regression. We are pretty used to linear regression, which requires at least 2 points and we connect them through a line. Like this:
To get anything between P1 and P2, we use this formula to interpolate the result¹. There’s polynomial regression, in which we try to pass through all points using a polynomial. Like a parabolic interpolation:
There are several ways to make regressions based on a discrete set of points, and we should select one to approximate our solution u. Preferably, we want something that has the piecewise polynomial form that we have “selected” for our approximation. Fortunately, there is such a regression model, which are Lagrange Polynomials.
The definition of Lagrange polynomials is: Given a set of points (xᵢ,yᵢ), with no values xᵢ equal, the Lagrange polynomial is the polynomial of lowest degree that assumes at each value xᵢ the corresponding value yᵢ. Taking it straight from Wikipedia, this means:
If you look closely, we managed to describe the function as a polynomial of x and a set of discrete points y, just like we did for our solution:
These lⱼ(x) polynomials are the basis of our approximation, and we will call them from now on the basis functions. So… how we calculate them? Where to sample a set of points in order to define the basis functions?
As the basis functions require points in space to be computed, the only reasonable place to look for a set of sampling points that is related to our problem is our domain Ω. So, we should discretize our domain as a set of points and make a set of Lagrange basis from this set of points, right?
Well, not exactly. There’s a big problem with polynomial regression of high degrees: oscillation. In order to accommodate the largest number of points inside this mathematical formula, when you get to high order polynomials, they start to behave like rollercoasters in order to correct the derivative and trajectory of the line. Again, from Wikipedia:
As you can notice from the Picture, the interpolation is close to perfection near the center. However, as we move from towards the edges, the interpolation greatly misses the original function. That’s the main problem with higher-order interpolations using Lagrange Polynomials.
So, how do we fix this issue? The answer is to keep the polynomial order low (p < 3) and work with small patches of the geometry. We then approximate the domain as a sum of these small patches²! We engineer really love approximating things. These small patches are what we call elements (big hints here).
Now, we have created a basis to describe our solution as a linear combination of. We still don’t know the coefficients that make up our solution (if we plug the physical values of the points, we will get the geometry, not the solution). As a linear combination whose only unknown are the coefficients, we should be able to create a system of equations to solve and find these coefficients — so we go back to the integral problem:
Now, let’s substitute our approximation of u as a sum of Lagrange polynomials:
The same assumptions that we have made for u will be made for v. We will describe it as a linear combination of our Lagrange polynomials with a different set of coefficients:
Thanks to the linear properties of the polynomials, we can isolate cᵢ from the problem, as they are constants. This will create a set of equations whose only unknowns are dⱼ:
We can then evaluate these integrals. Several ways of evaluating integrals exist (analytical, numerical, quadrature, etc). Most commercial FEM software use either a pre-processed formula (as they know the polynomials and their derivatives beforehand) or using a quadrature rule (which is exact for polynomials). Either way, as one can deduct by the subscripts, you will end with a matrix system in the form of:
Which can be solved with a matrix inversion of K. Historically, we call K the stiffness matrix, F the Force vector and d the displacement. This is due to the origins of the method in linear elasticity problems. After solving for dⱼ, we plug them back into the polynomial approximation to find our u.
And there we go, people. We have just described how to do a numerical simulation. There are still some very important issues with this approach that we haven’t mentioned, which will be the topic of next week: boundary conditions. Stay tuned! :)
¹ In general, for a statistical linear regression model, you want more points than just 2. You then create the line that minimizes the distance to each point of the set. You will (almost) never have a perfect line connecting all points, but you will have the best line fit.
²In order for this sum of small patches to work, we still have to define the Lagrange polynomials to be equal to zero outside the patch, otherwise, we wouldn’t get the domain as a sum of all patches.
