Solvers and Time Integration
Various implementations of the method of lines to progress through time. The goal is to implement the code in python and not rely on existing solvers.
Bram Lagerweij COHMAS Mechanical Engineering KAUST 2021
- solvers.backwardEuler(pde, u, dt, t_end)
Itterate a through time with the backward Eurler method.
Lets assume that, through any type of discretization, the time derivative was obtained. This time derivative can be represented through linear algabra as:
\[M\,u_t = K\,u + b \qquad \text{that is} \qquad u_t = M^{-1}(K\,u + b)\]where \(M\) is the mass matrix, \(K\) the siffness and transport matrix and vector \(b\) the right hand side. these are obtained from approximations of the spatial derivatives defined by the functien provided to func
The backward Euler method predicts the field of our function based upon information of the previous timestep only. Imagine that we are at timestep \(n\) and want to predict our field at timestep \(u^{(n+1)}\). Now a backward finite difference approximation used the time derivative of the next timestep, wich is not yet known:
\[u^{(n+1)}_t = \frac{ -u^{(n)} + u^{(n+1)} }{dt}\]That is we can predict our field in the future timestep as:
\[u^{(n+1)} = u^{(n)} + dt\, u^{(n+1)}_t\]in which we substitute the linear algabra representation of our PD.
\[u^{(n+1)} = u^{(n)} + dt\, M^{-1}(K u^{n+1} + b)\]It is important to notic that there is a term with an unknown, as that is at time step :math:`n+1’ on both sides of the equation. Now we rewrite it into a system of equations where we find all unknowns on the left hand side and all knownn on the right hand side.
\[(M - dt\,K)\,u^{(n+1)} = M\,u^{(n)} + dt b\]This is a system of equations which can be solved.
- Parameters
pde (tuple) – The linear algabra objects of the pde \(M\,u_t = K\,u + b\).
u (array_like) – The field at the start \(u(t=0)\).
dt (float) – The size of the time step.
t_end (float) – Time at termination.
- Returns
The function for all time steps.
- Return type
array_like
- solvers.forwardEuler(pde, u, dt, t_end)
Itterate a through time with the forward Eurler method.
Lets assume that, through any type of discretization, the time derivative was obtained. This time derivative can be represented through linear algabra as:
\[M \, u_t = K \, u + b \qquad \text{that is} \qquad u_t = M^{-1}(K u + b)\]where \(M\) is the mass matrix, \(K\) the siffness and transport matrix and vector \(b\) the right hand side. these are obtained from approximations of the spatial derivatives defined by the functien provided to func.
The backward Euler method predicts the field of our function based upon information of the previous timestep only. Imagine that we are at timestep \(n\) and want to predict our field at timestep \(u^{(n+1)}\). Now a forward finite difference approximation is used:
\[u^{(n)}_t = \frac{-u^{(n)} + u^{(n+1)} }{dt}\]That is we can predict our field in the future timestep as:
\[u^{(n+1)} = u^{(n)} + dt\, u^{(n)}_t\]Now from our linear algabra implementation we substitute \(u_t\)
\[u^{(n+1)} = u^{(n)} + dt\, M^{-1}(K u^{(n)} + b)\]most important of all is to see that everything on the right hand side is exactly known. Thus the updated field can be calculated directly. However For this purpouse we would have to invert the mass matrix. If the mass matrix is the identity matrix this is simple, but in generally this is not the case. As we don’t want to invert large matrices, we multiply all terms by \(M\).
\[M u^{(n+1)} = M u^{(n)} + dt\,(K u^{(n)} + b)\]Which is a system of equations as everything on the right hand side is known and can be calculated directly.
Notes
This code will recognize if \(M\) is the identity matrix and, in that case it will solve the problem directly, avoiding the need to solve a sytem of equations.
- Parameters
pde (tuple) – The linear algabra objects of the pde \(M\,u_t = K\,u + b\).
u (array_like) – The field at the start \(u(t=0)\).
dt (float) – The size of the step.
t_end (float) – Time at termination.
- Returns
The function for all time steps.
- Return type
array_like
- solvers.solve(K, b)
Solve a time independed problem.
- Parameters
func (callable) – The linear algabra problem that we want to solve \(K\,u = b\).
args (tuple, optional) – The parameters into the PDE approximation. Defealts to an empty tuple.
- Returns
The vector containing \(u\).
- Return type
array_like