# A Compact FEM Implementation for Parabolic Integro-Differential Equations in 2D
## 7 Overview of an Implementation in Matlab  
A minimal MATLAB implementation of the presented algorithm is available from Supplementary Materials. There are six zipfiles which contain the necessary MATLAB files therein:  
* BE_LRR_unvectorized.zip (backward-Euler, left rectangular rule, unvectorized);
* BE_LRR_vectorized.zip (backward-Euler, left rectangular rule, vectorized);
* BE_RRR_unvectorized.zip (backward-Euler, right rectangular rule, unvectorized);
* BE_RRR_vectorized.zip (backward-Euler, right rectangular rule, vectorized);
* CN_unvectorized.zip (Crank-Nicolson, unvectorized);
* CN_vectorized.zip (Crank-Nicolson, vectorized). ```

**Algorithm 11** Assembling the mass matrix with vectorization. To run the demo program, uncompress the desired zip file and run the Main.m file; this:* Specifies the geometry;
* Computes the solution;
* Plots the solution. In particular, the computation of the solution is carried out by SCHEME_NAME_PIDE.m, which sets the initial conditions via U0.m, and assembles the stiffness matrix and the mass matrix via StiffAssemb.m and MassAssemb.m, respectively; thereafter, solve_SCHEME_NAME.m updates the boundary conditions via Ud.m, assembles the load vector via f.m and obtains the solution at the current time step. Furthermore, note that whilst most of the subroutines are the same in both the vectorized and unvectorized versions of the code, those for StiffAssemb.m and MassAssemb.m are necessarily slightly different, in line with the discussion in Sections 5 and 6. ## 8 Numerical Results

In this section, we study a numerical example and compare the results obtained using different versions of the MATLAB-based sets of code; in addition, solutions were obtained using COMSOL Multiphysics. For the example, we consider the L-shaped domain \(\Omega=[-1,1]^{2}\backslash[-1,0]\times[0,1]\) mentioned earlier; take

\[f(x,y,t)=(1-2t)\pi^{2}u(x,y,t) \tag{25}\]

and

\[B(t,s)=\exp(-\pi^{2}(t-s)), \tag{26}\]

which is a kernel that occurs in different settings in many of the references given earlier, and others [2; 3; 4; 7; 14; 53; 54]; this leads to the following exact solution for \(u\):

\[u(x,y,t)=\exp(-\pi^{2}t)\sin(\pi x)\sin(\pi y). \tag{27}\]

Figure 2a shows the exact solution to the problem, whereas Figure 2b corresponds to the backward-Euler approximation, where the left rectangular rule is applied to treat the Volterra integral term. Figure 2c,d show the solutions obtained using the backward-Euler approximation with the right rectangular rule and the Crank-Nicolson approximation with the trapezoidal rule, respectively, to treat the Volterra integral term. In addition, Figure 3 shows a comparison of the analytical solution with the solution computed by COMSOL Multiphysics with 66142 P1 elements for \(t=T,y=-0.5\), where \(T=0.1\); as expected, the agreement is very good. To illustrate the advantage of the vectorized code for each numerical scheme over the unvectorized one, which uses _for_ loops for the assembly process, we have presented a comparison of the runtimes with respect to the number of elements, \(N\). This is given in Table 3, which also shows the runtimes for COMSOL Multiphysics. Note that all sets of code were run on an Intel Core i7 Notebook with a 1.7 GHz processor and 8 GB of RAM. We can note from Table 3 that vectorized versions of the MATLAB sets of code presented in this article on the finest mesh get close to COMSOL's performance. Nevertheless, there are several factors behind the results. On the one hand, in the developed MATLAB sets of code, we discretize the Volterra integral term; as a result, one has to store the memory term values at all of the time discretization points. On the other hand, because of the particular form of the kernel in (26), and as explained in Appendix C, COMSOL Multiphysics does not need to evaluate the integral at all; even so, the difference in computational time is not decisive.