# A Compact FEM Implementation for Parabolic Integro-Differential Equations in 2D
## 7 Overview of an Implementation in Matlab
* **Code for the backward-Euler scheme with the right rectangular rule, corresponding to Equation (9) and employing (8):**

```
1:procedureBERRFEMPIDE(\(c4n,n4e,dir,N,dt\))
2:\(r\leftarrow\)NROW(\(c4n\))
3:\(U\leftarrow\)ZEROS(\(r,N+1\))
4:\(FN\leftarrow\)UNIQUE(\([1\ 2\ \ldots\ r])\backslash\)UNIQUE(\(dir\))
5:\(A\leftarrow\)StiffAssemb(\(n4e,c4n\))
6:\(B\leftarrow\)MassAssemb(\(n4e,c4n\))
7:\(U(:,1)\leftarrow\)U0(\(c4n\))
8:for\(n\gets 2\) to \((N+1)\),do
9:\(U(:,n)\leftarrow\)solveBERRR(\(c4n,n4e,dir,dt,FN,A,B.n,U\))
10:endfor
11:return\(U\)
12:endprocedure
```

**Algorithm A2** Backward-Euler scheme with right rectangular rule. ```
1:proceduresolveBERRR(\(c4n,n4e,dir,dt,FN,A,B,n,U\))
2:\(r\leftarrow\)NROW(\(c4n\))
3:\(b\leftarrow\)SPARSE(\(r,1\))
4:for\(j\gets 1\) to NROW(\(n4e\)),do
5:\(M_{1}\leftarrow[[1\ 1\ 1];c4n(n4e(j,:),:)^{T}]\)
6:\(S\leftarrow\)SUM(\(c4n(n4e(j,:),:)\))
7:\(b(n4e(j,:))\gets b(n4e(j,:))+\)DET(\(M_{1}\))\(*dt*\frac{1}{6}*\)F(\(\frac{5}{3},(n-1)*dt\))
8:endfor
9:\(Int\gets 0\)
10:if\(n>2\),then
11:\(Int\gets dt^{2}*A*U(:,1)*\textsc{intgrad}((n-1)*dt,dt)\)
12:for\(k\gets 2\) to \((n-2)\),do
13:\(Int\gets Int+\big{(}dt^{2}*A*U(:,k)*\textsc{intgrad}((n-1)*dt,k*dt)\big{)}\)
14:endfor
15:endif
16:\(b\gets b+(B*U(:,n-1))+Int\)
17:\(u\leftarrow\)SPARSE(\(r,1\))
18:\(u_{d}\leftarrow\)UNIQUE(\(dir\))
19:\(u(u_{d})\leftarrow\)UD(\(c4n(u_{d},:),(n-1)*dt\))
20:\(b\gets b-(dt*A+B+dt^{2}*A*\textsc{intgrad}((n-1)*dt,(n-1)*dt))*u\)
21:\(M_{1}\gets dt*A(FN,FN)+B(FN,FN)+dt^{2}*A(FN,FN)*\textsc{intgrad}((n-1)*dt,(n-1)*dt)\)
22:\(u(FN)\leftarrow\)INV(\(M_{1}\))\(*b(FN\))
23:\(solBERR\gets u\)
24:return\(solBERR\)
25:endprocedure
```

**Algorithm A3** Backward-Euler scheme with right rectangular rule. * **Code for Crank-Nicolson scheme with trapezoidal rule, corresponding to Equation (14):**

```
1:procedureCNFEMPIDE(\(c4n,n4e,dir,N,dt\))
2:\(r\leftarrow\)NROW(\(c4n\))
3:\(U\leftarrow\)ZEROS(\(r,N+1\))
4:\(FN\leftarrow\)UNIQUE(\([12\ \dots\ r])\backslash\)UNIQUE(\(dir\))
5:\(A\leftarrow\)StiffAssemb(\(n4e,c4n\))
6:\(B\leftarrow\)MassAssemb(\(n4e,c4n\))
7:\(U(:,1)\leftarrow\)U0(\(c4n\))
8:for\(j\gets 2\) to \((N+1)\),do
9:\(U(:,j)\leftarrow\)solveCN(\(c4n,n4e,dir,dt,FN,n,U,A,B\))
10:endfor
11:return\(U\)
12:endprocedure
```

**Auxiliary function**

```
1:proceduresolveCN(\(c4n,n4e,dir,dt,FN,n,U,A,B\))
2:\(r\leftarrow\)NROW(\(c4n\))
3:\(b\leftarrow\)SPARSE(\(r,1\))
4:for\(j\gets 1\) to NROW(\(n4e\)),do
5:\(X_{1}\leftarrow\)DET(\([11\ 1];c4n(n4e(j,:),:)^{T}]\))
6:\(X_{2}\leftarrow\)SUM(\(c4n(n4e(j,:),:)\))
7:\(X_{3}\leftarrow\)SUM(\(c4n(n4e(j,:),:)\))
8:\(b(n4e(j,:))\gets b(n4e(j,:))+X_{1}*dt*\frac{1}{12}*(\textsc{f}(X_{2}/3,(n-1)* dt)+\textsc{f}(X_{3}/3,(n-2)*dt))\)
9:endfor
10:\(Int\leftarrow\frac{1}{8}. *dt^{2}*A*U(:,1)*(2*\textsc{integral}((n-1.5)*dt,0)+ \textsc{integral}((n-1.5)*dt,(n-1.5)*dt))\)
11:if\(n>2\),then
12:\(Int\leftarrow\frac{1}{2}*dt^{2}*A*U(:,1)*\textsc{integral}((n-1.5)*dt,0)\)\(+\frac{1}{8}*dt^{2}*A*U(:,n-1)*(\textsc{integral}((n-1.5)*dt,(n-1.5)*dt)\)\(+6*\textsc{integral}((n-1.5)*dt,(n-2)*dt))\)
13:for\(k\gets 2\) to \((n-2)\),do
14:\(Int\gets Int+dt^{2}*A*U(:,k)*\textsc{integral}((n-1.5)*dt,(k-1)*dt)\)
15:endfor
16:endif
17:\(b\gets b+(B-\frac{dt}{2}*A)*U(:,n-1)+Int\)
18:\(u\leftarrow\)SPARSE(\(r,1\))
19:\(u_{d}\leftarrow\)UNIQUE(\(dir\))
20:\(u(u_{d})\leftarrow\)UD(\(c4n(u_{d},:),(n-1)*dt\))
21:\(b\gets b-\left(\left(\frac{dt}{2}-\frac{dt^{2}}{8}*\textsc{integral}((n-1.5)*dt,(n-1)*dt) \right)*A+B\right)*u\)
22:\(M_{1}\leftarrow\left(\left(\frac{dt}{2}-\frac{dt^{2}}{8}*\textsc{integral}((n-1.5)*dt,(n-1)*dt) \right)*A(FN,FN)+B(FN,FN)\right)\)
23:\(u(FN)\leftarrow\)INV(\(M1)*b(FN\))
24:\(solCN\gets u\)
25:return\(solCN\)
26:endprocedure
```

**Algorithm A3** Crank-Nicolson scheme with the trapezoidal rule. ## Appendix B Vectorized Pseudocode

In this section, the vectorized sets of pseudocode are presented.