微分方程数值解 第六章（3）

二维问题的矩形元

一、问题与变分形式

二维 Poisson 方程（ Dirichlet 边界）：
$$
\begin{cases}
-\Delta u = f(x,y), & (x,y) \in \Omega = [a,b] \times [c,d] \
u|_{\partial\Omega} = g(x,y)
\end{cases}
$$
变分形式（虚功原理）：

找 $u \in H_0^1(\Omega)$（满足边界条件的 Sobolev 空间 ），使
$$
\int_{\Omega} \nabla u \cdot \nabla v , dxdy = \int_{\Omega} fv , dxdy
$$
对任意 $v \in H_0^1(\Omega)$ 成立（$\Delta$ 是 Laplace 算子，$\nabla$ 是梯度 ）

二、Galerkin 离散（有限元方法）

有限元空间 $U_h$：

取 $U_h$ 为分片双线性函数空间（Lagrange 型，在矩形剖分单元上线性 ），满足 $u_h|_{\partial\Omega} = g_h$（$g_h$ 是 $g$ 的边界插值 ）

离散方程：

找 $u_h \in U_h$ 使
$$
\int_{\Omega} \nabla u_h \cdot \nabla v_h , dxdy = \int_{\Omega} f v_h , dxdy
$$
对任意 $v_h \in U_h$（$v_h|_{\partial\Omega} = 0$ ）成立，代入基函数展开得线性方程组

三、矩形剖分与基函数

区域剖分：

将 $\Omega$ 剖分为矩形单元 $R_{ij} = [x_i, x_{i+1}] \times [y_j, y_{j+1}]$ ，步长 $\Delta x = x_{i+1}-x_i$，$\Delta y = y_{j+1}-y_j$ 。

双线性基函数：

单元上基函数为 $1, x, y, xy$ 的线性组合（双线性插值 ），通过节点值 $u_{ij}, u_{i+1,j}, u_{i,j+1}, u_{i+1,j+1}$ 构造，满足分片双线性与 Kronecker delta 性质（节点处取 1，其余为 0 ）

标准化与坐标变换

对矩形单元 $R_{ij} = [x_i, x_{i+\Delta x}] \times [y_j, y_{j+\Delta y}]$，通过仿射变换标准化到参考单元 $[0,1] \times [0,1]$：
$$
\xi = \frac{x - x_i}{\Delta x}, \quad \eta = \frac{y - y_j}{\Delta y}
$$

双线性插值函数构造

设单元顶点值为 $u_{00}, u_{10}, u_{01}, u_{11}$（对应 $(\xi,\eta)=(0,0),(1,0),(0,1),(1,1)$ ），双线性插值函数 $p(\xi,\eta)$ 满足：
$$
p(\xi,\eta) = c_0 + c_1 \xi + c_2 \eta + c_3 \xi\eta
$$
通过顶点插值条件
$$
p(0,0)=u_{00},\quad p(1,0)=u_{10}\
p(0,1)=u_{01},\quad p(1,1)=u_{11}
$$
解得 Lagrange 插值多项式：
$$
p(\xi, \eta) = u_{00}(1 - \xi)(1 - \eta) + u_{10} \xi(1 - \eta) + u_{01}(1 - \xi)\eta + u_{11} \xi \eta
$$

二维 Lagrange 基函数：
$$
\begin{align*}
p^{(0,0)}(\xi, \eta) &= (1 - \xi)(1 - \eta) = N_0(\xi)N_0(\eta) \
p^{(1,0)}(\xi, \eta) &= \xi(1 - \eta) = N_1(\xi)N_0(\eta) \
p^{(0,1)}(\xi, \eta) &= (1 - \xi)\eta = N_0(\xi)N_1(\eta) \
p^{(1,1)}(\xi, \eta) &= \xi \eta = N_1(\xi)N_1(\eta)
\end{align*}
$$
表示为Lagrange 基函数的乘积形式：
$$
p(\xi,\eta) = u_{00} N_0(\xi) N_0(\eta) + u_{10} N_1(\xi) N_0(\eta) + u_{01} N_0(\xi) N_1(\eta) + u_{11} N_1(\xi) N_1(\eta)
$$
基函数为：
$$
N_0(t) = 1 - t, \quad N_1(t) = t
$$

四、单元插值函数（$ (x,y) \in R_{ij} $ ）

$$
u_h(x,y) = u_{ij} p_{ij}^{(0,0)}(\xi,\eta) + u_{i+1,j} p_{ij}^{(1,0)}(\xi,\eta) + u_{i,j+1} p_{ij}^{(0,1)}(\xi,\eta) + u_{i+1,j+1} p_{ij}^{(1,1)}(\xi,\eta)
$$

给矩形单元排序（从左下至上依次记为 $R_{11},R_{12},\cdots,R_{MN}$）
$$
u_h(x,y) =
\begin{cases}
u_{11} p_{11}^{(0,0)}(\xi,\eta) + u_{21} p_{11}^{(1,0)}(\xi,\eta) + u_{12} p_{11}^{(0,1)}(\xi,\eta) + u_{22} p_{11}^{(1,1)}(\xi,\eta), \quad R_{11}\
u_{12} p_{12}^{(0,0)}(\xi,\eta) + u_{22} p_{12}^{(1,0)}(\xi,\eta) + u_{13} p_{12}^{(0,1)}(\xi,\eta) + u_{23} p_{11}^{(1,1)}(\xi,\eta), \quad R_{12} \
\quad \quad \vdots \
u_{1,N} p_{1N}^{(0,0)}(\xi,\eta) + u_{2,N} p_{1N}^{(1,0)}(\xi,\eta) + u_{1,N+1} p_{1N}^{(0,1)}(\xi,\eta) + u_{2,N+1} p_{1N}^{(1,1)}(\xi,\eta) , \quad R_{1N}\
\quad \quad \vdots \
u_{M,N} p_{MN}^{(0,0)}(\xi,\eta) + u_{M+1,N} p_{MN}^{(1,0)}(\xi,\eta) + u_{M,N+1} p_{MN}^{(0,1)}(\xi,\eta) + u_{M+1,N+1} p_{MN}^{(1,1)}(\xi,\eta), \quad R_{MN}
\end{cases}
$$

五、整体双线性基函数

$$
\varphi_{ij}(x,y) =
\begin{cases}
\left(1 - \frac{|x - x_i|}{\Delta x}\right)\left(1 - \frac{|y - y_j|}{\Delta y}\right), & (x,y) \in R_{ij} \cup R_{i+1,j} \cup R_{i,j+1} \cup R_{i+1,j+1} \
0, & \text{其他}
\end{cases}
$$

六、自由度与空间维数

若把 $u_h$ 在所有节点的函数值都当作自变量，节点数 $ (M+1) \times (N+1) $，故：
$$
\text{$\varphi_{i,j}$ 的个数为} (M+1)(N+1), \quad \dim(U_h) = (M+1)(N+1)
$$

七、网格剖分与数据结构

一、网格剖分与标号规则

节点标号：

矩形区域剖分为 $M \times N$ 单元，节点总数 $K = (M+1)(N+1)$，按顺序标号 $1, 2, \dots, K$ 。

单元标号：
单元总数 $M \cdot N$，按顺序标号 $1, 2, \dots, M \cdot N$ ，每个单元对应 4 个顶点（左下→右下→右上→左上 ）。

二、数据结构定义

节点数组（node）：

存储节点坐标，维度 $K \times 2$ ，每行对应一个节点的 $(x,y)$ 坐标：
$$
\text{node} = \begin{bmatrix} x^{(1)} & y^{(1)} \ \vdots & \vdots \ x^{(K)} & y^{(K)} \end{bmatrix}
$$
单元数组（elem）：
存储单元顶点的节点标号，维度 $M \cdot N \times 4$ ，每行对应一个单元的 4 个顶点标号（按左下→右下→右上→左上顺序 ）：
$$
\text{elem} = \begin{bmatrix} i_{j,1} & i_{j,2} & i_{j,3} & i_{j,4} \ \vdots & \vdots & \vdots & \vdots \end{bmatrix}
$$

$$
R_{ij}(i,j)=[u_{ij},u_{i+1,j},u_{i,j+1},u_{i+1,j+1}]
$$

八、单元刚度矩阵积分

$$
\iint_{R_{ij}} \nabla u_h \cdot \nabla v_h , dxdy = \iint_{R_{ij}} \left( \frac{\partial u_h}{\partial x} \frac{\partial v_h}{\partial x} + \frac{\partial u_h}{\partial y} \frac{\partial v_h}{\partial y} \right) dxdy
$$

九、插值函数与导数（局部坐标 $ \xi, \eta $ ）

单元内 $ u_h(x,y) $ 展开：
$$
u_h(x,y) = u_{ij} N_0(\xi)N_0(\eta) + u_{i+1,j} N_1(\xi)N_0(\eta) + u_{i,j+1} N_0(\xi)N_1(\eta) + u_{i+1,j+1} N_1(\xi)N_1(\eta)
$$
$ x $ 方向导数（$ \displaystyle \frac{\partial \xi}{\partial x} = \frac{1}{\Delta x} $ ）：
$$
\frac{\partial u_h}{\partial x} = \frac{1}{\Delta x} \left[ u_{ij} (- (1 - \eta)) + u_{i+1,j} (1 - \eta) + u_{i,j+1} (- \eta) + u_{i+1,j+1} \eta \right]
$$
$v_h$ 同理

十、单元刚度矩阵计算（$ x $ 方向示例）

$$
\begin{align}
\iint_{R_{ij}} \frac{\partial u_h}{\partial x} \frac{\partial v_h}{\partial x},dxdy&= \int_{R_{ij}} [u_{ij}[-(1-\eta)]+\cdots]\cdot[v_{ij}[-(1-\eta)]+\cdots]\cdot\frac{1}{\Delta x_i^2},dxdy\
&=\int_0^1 \int_0^1 [u_{ij}[-(1-\eta)]+\cdots]\cdot[v_{ij}[-(1-\eta)]+\cdots]\cdot\frac{1}{\Delta x_i^2}\cdot \Delta x\cdot\Delta y,d\xi d\eta\
&=(v_{ij},v_{i+1,j}, v_{i,j+1},v_{i+1,j+1}) \begin{bmatrix}
\displaystyle \int_0^1 (1-\eta)^2 d\eta &\displaystyle \int_0^1 -(1-\eta)^2 d\eta & \times &\times \
\times & \times & \times &\times \
\times & \times & \times &\times \
\times & \times & \times &\times
\end{bmatrix} \frac{\Delta y}{\Delta x} \begin{bmatrix} u_{i j} \ u_{i+1,j} \ u_{i,j+1}\u_{i+1,j+1} \end{bmatrix}\
&=(v_{ij},\cdots,v_{i+1,j+1})K_ex\begin{bmatrix} u_{i j} \ \vdots\u_{i+1,j+1} \end{bmatrix}
\end{align}
$$

$K_ey$ 同理

十一、整体单元刚度矩阵

$$
K_e = K_{e}x + K_{e}y
$$

$ K_{e}y $ 为 $ y $ 方向单元刚度矩阵，对称推导

十二、总刚度矩阵 $K$ 的组装

初始化：

总刚度矩阵 $K$ 维度为 $(M+1)(N+1) \times (M+1)(N+1)$，初始化为零矩阵：
$$
K = \text{zeros}(((M+1)(N+1), (M+1)(N+1)))
$$
单元循环组装：

遍历每个单元（标号 $k=1$ 到 $ne = M \cdot N$ ），计算单元刚度矩阵并叠加到总矩阵对应位置（通过节点标号索引 ）

伪代码示意 ：
$$
\text{for } k = 1 : n_e \
\quad \text{将单元刚度矩阵按节点索引叠加到 } K \
\text{end}
$$

十三、载荷向量的积分计算

对载荷项 $\displaystyle \int_{\Omega} f v_h , dxdy$，在矩形单元 $R_{ij}$ 上通过双线性基函数展开，利用变量替换（标准化到参考单元 $[0,1] \times [0,1]$ ），将积分转化为：
$$
\begin{align}
\int_{R_{ij}} f v_h , dxdy &=\int_0^1\int_0^1f(x_i+\Delta x\cdot\xi,y_j+\Delta y\cdot\eta)[v_{ij}N_0(\xi)N_0(\eta)+\cdots+v_{i+1,j+1}N_1(\xi)N_1(\eta)]\cdot \Delta x\cdot\Delta y,d\xi d\eta\

&= (v_{ij}, \dots, v_{i+1,j+1}) \begin{bmatrix} \displaystyle \int_0^1 \int_0^1 f N_0(\xi)N_0(\eta) d\xi d\eta \ \vdots \ \displaystyle\int_0^1 \int_0^1 f N_1(\xi)N_1(\eta) d\xi d\eta \end{bmatrix} \Delta x \Delta y\
&\approx(v_{ij}, \dots, v_{i+1,j+1})\begin{bmatrix} \frac{1}{4}f_{ij} \ \vdots \\frac{1}{4}f_{i+1,j+1}\end{bmatrix}\Delta x \Delta y
\end{align}
$$
$\Delta x, \Delta y$ 是单元步长，$v_j$ 是节点载荷系数

十四、最终变分形式

离散后的变分方程表示为：
$$
\mathbf{V}^T A \mathbf{U} = \mathbf{V}^T \mathbf{F}\
A=K
$$
其中，$A = K$ 是总刚度矩阵，$\mathbf{U}$ 是节点值向量，$\mathbf{V}$ 是检验函数向量，$\mathbf{F}$ 是载荷向量，体现有限元离散的线性方程组形式（$K \mathbf{U} = \mathbf{F}$ ）

十五、边界条件处理

节点分类：

边界节点 $ \text{bdnode} $（已知 $ U(\text{bdnode}) $ ）、自由节点 $ \text{freenode} $（待求 $ U(\text{freenode}) $ ）

变分约束：

检验函数 $ V(\text{bdnode}) = 0 $，自由节点 $ V(\text{freenode}) $ 任意，由此从整体系统 $ A U = F $ 中提取自由节点子系统：
$$
A(\text{freenode}, :) U - F(\text{freenode}) = 0
$$
子系统分解：
展开为自由节点与边界节点的矩阵块形式：
$$
A(\text{freenode}, \text{freenode}) U(\text{freenode}) + A(\text{freenode}, \text{bdnode}) U(\text{bdnode}) = F(\text{freenode})
$$
求解自由节点：
移项得自由节点的代数方程：
$$
A(\text{freenode}, \text{freenode}) U(\text{freenode}) = F(\text{freenode}) - A(\text{freenode}, \text{bdnode}) U(\text{bdnode})
$$
最终解为：
$$
U(\text{freenode}) = A(\text{freenode}, \text{freenode})^{-1} \left( F(\text{freenode}) - A(\text{freenode}, \text{bdnode}) U(\text{bdnode}) \right)
$$
核心是有限元中Dirichlet 边界条件的离散处理方法，通过节点分类和矩阵分块，将已知边界值代入整体系统，求解自由节点未知量

三角元法

一、问题与变分形式

二维 Poisson 方程（ Dirichlet 边界）：
$$
\begin{cases}
-\Delta u = f(x,y), & (x,y) \in \Omega \
u|_{\partial\Omega} = g(x,y)
\end{cases}
$$
变分形式（虚功原理）：

找 $u \in H_0^1(\Omega)$（满足边界条件的 Sobolev 空间 ），使
$$
\int_{\Omega} \nabla u \cdot \nabla v , dxdy = \int_{\Omega} fv , dxdy
$$
对任意 $v \in H_0^1(\Omega)$ 成立（$\Delta$ 是 Laplace 算子，$\nabla$ 是梯度 ）

二、三角元法的核心优势

三角形剖分可逼近任意复杂曲线边界（通过多边形域近似 $\Omega$ ），适配工程中不规则区域（如飞机翼型、地质构造 ）的数值求解

三、三角剖分规则

将多边形域 $\Omega$ 剖分为有限个三角形单元 $\Delta_j$，需满足：

1. 单元不重叠（$\Delta_i^o \cap \Delta_j^o = \emptyset, i \neq j$ ，$^o$ 表示内部 ）

2. 无悬挂节点（三角形顶点不落在另一三角形的边上 ）

进一步定义正则剖分（最小角有下界 ）和拟一致剖分（单元尺寸均匀 ），保证有限元解的收敛性

四、Lagrange 型三角元空间

有限元空间 $U_h$ 定义为：
$$
U_h = \left{ u \in H_0^1(\Omega) \mid u|{\Delta_j} \in P^m(\Delta_j) \right}
$$
其中，$P^m(\Delta_j)$ 是三角形单元上的 $m$ 次 Lagrange 多项式，形式为：
$$
P_m(x,y) = \sum{i+j \leq m} c_{ij} x^i y^j
$$
$m$ 次多项式的系数个数为 $\frac{1}{2}(m+1)(m+2)$，需取对应数量的插值节点保证唯一性

五、单元与插值函数推导

对二维 Poisson 方程的三角剖分单元 $\Delta_{(123)}$（节点按逆时针排序为 $(x_1,y_1),(x_2,y_2),(x_3,y_3)$ ），构造线性插值函数 $p(x,y) = ax + by + c$ ，满足顶点插值条件 $p(x_i,y_i) = u_i$（$i=1,2,3$ ）

通过顶点条件建立线性方程组：
$$
\begin{cases}
a x_1 + b y_1 + c = u_1 \
a x_2 + b y_2 + c = u_2 \
a x_3 + b y_3 + c = u_3
\end{cases}
$$
写成矩阵形式 $A \begin{bmatrix} a \ b \ c \end{bmatrix} = \begin{bmatrix} u_1 \ u_2 \ u_3 \end{bmatrix}$ ，其中 $A = \begin{bmatrix} x_1 & y_1 & 1 \ x_2 & y_2 & 1 \ x_3 & y_3 & 1 \end{bmatrix}$

解得系数：
$$
\begin{bmatrix} a \ b \ c \end{bmatrix} = A^{-1} \begin{bmatrix} u_1 \ u_2 \ u_3 \end{bmatrix}
$$
因此，插值函数可表示为：
$$
p(x,y) = (x, y, 1) A^{-1} \begin{bmatrix} u_1 \ u_2 \ u_3 \end{bmatrix}
$$

六、面积坐标（Lagrange 基函数）

定义面积坐标 $L_1,L_2,L_3$：
$$
(L_1, L_2, L_3) = (x, y, 1) A^{-1}
$$
则插值函数简化为：
$$
p(x,y) = u_1 L_1(x,y) + u_2 L_2(x,y) + u_3 L_3(x,y)
$$
其中，$L_i(x,y)$ 是线性 Lagrange 基函数，满足 $L_i(x_j,y_j) = \delta_{ij}$（Kronecker delta ），且 $L_1 + L_2 + L_3 = 1$
$$
\begin{align}
L_1&=\frac{1}{\det(A^T)}\begin{vmatrix} x &x_2 &x_3 \ y & y_2 & y_3 \ 1 & 1 & 1 \end{vmatrix}\
L_2&=\frac{1}{\det(A^T)}\begin{vmatrix} x_1 &x &x_3 \ y_1 & y & y_3 \ 1 & 1 & 1 \end{vmatrix}\
L_3&=\frac{1}{\det(A^T)}\begin{vmatrix} x_1 &x_2 &x \ y_1 & y_2 & y \ 1 & 1 & 1 \end{vmatrix}
\end{align}
$$

七、几何意义（Cramer 法则）

通过 Cramer 法则计算 $L_i$，发现分母 $ \det(A^T) = \det(A) $ 与三角形面积相关：
$$
\det(A) = 2 S_{\Delta_{(123)}}
$$
其中
$$
S_{\Delta_{(123)}} = \frac{1}{2} \left| \det \begin{bmatrix} x_1 & y_1 & 1 \ x_2 & y_2 & 1 \ x_3 & y_3 & 1 \end{bmatrix} \right|
$$
记 $(x_1,y_1)\rightarrow(x_2,y_2)$ 的向量为 $\vec{a}$，$(x_1,y_1)\rightarrow(x_3,y_3)$ 的向量为 $\vec{b}$，则
$$
\vec a \times \vec b=\begin{vmatrix} i &j & k \ x_2-x_1 & y_2-y_1 & 0 \ x_3-x_1 & y_3-y_1 & 0 \end{vmatrix}
$$

$$
S_{\Delta_{(123)}} = \frac{1}{2} \begin{vmatrix}\vec a \times \vec b\end{vmatrix}
$$

八、面积坐标与直角坐标的关系

$$
\begin{align}
L_1&=\frac{2S_1}{2S}=\frac{S_1}{S}=\frac{(y_2 - y_3)(x - x_3) + (x_3 - x_2)(y - y_3)}{(y_2 - y_3)(x_1 - x_3) + (x_3 - x_2)(y_1 - y_3)}\
L_2&=\frac{2S_2}{2S}=\frac{S_2}{S}=\frac{(y_3 - y_1)(x - x_3) + (x_1 - x_3)(y - y_3)}{(y_2 - y_3)(x_1 - x_3) + (x_3 - x_2)(y_1 - y_3)}\
L_3&=\frac{2S_3}{2S}=\frac{S_3}{S}=\frac{(y_1 - y_2)(x - x_2) + (x_2 - x_1)(y - y_2)}{(y_2 - y_3)(x_1 - x_3) + (x_3 - x_2)(y_1 - y_3)}\
\end{align}
$$

$$
x = L_1 x_1 + L_2 x_2 + L_3 x_3\
y = L_1 y_1 + L_2 y_2 + L_3 y_3
$$

九、单元插值函数

三角形单元上的线性插值函数表示为：
$$
u_h(x,y) = u_1 L_1(x,y) + u_2 L_2(x,y) + u_3 L_3(x,y)
$$
$u_i$ 是顶点 $i$ 的函数值，$L_i$ 是面积坐标基函数

十、全域基函数与自由度

通过 “单元拼接” 构造全域基函数 $\varphi_j(x,y)$（在节点 $P_j$ 处取 1，其他单元为 0 ），有限元空间 $U_h$ 的维度为节点总数 $nd$（自由度 ）
$$
u_h(x,y)=\sum_{j=1}^{nd}u_j\cdot \varphi_j(x,y)
$$

十一、单元刚度矩阵

单元刚度矩阵 $K_e$ 由变分形式的积分推导：
$$
\begin{align}
\int_{\Delta_j} \nabla u_h \cdot \nabla v_h , dxdy &= \int_{\Delta_j}(u_1^{(j)}\nabla L_1^{(j)}+u_2^{(j)}\nabla L_2^{(j)}+u_3^{(j)}\nabla L_3^{(j)})(v_1^{(j)}\nabla L_1^{(j)}+v_2^{(j)}\nabla L_2^{(j)}+v_3^{(j)}\nabla L_3^{(j)})\
&=(v_1^{(j)},v_2^{(j)}, v_3^{(j)}) \begin{bmatrix}
\displaystyle \int_{\Delta j} \nabla L_1^{(j)}\nabla L_1^{(j)} dxdy &\displaystyle \int_{\Delta j} \nabla L_1^{(j)}\nabla L_2^{(j)} dxdy & \times &\
\times & \times & \times \
\times & \times & \times \
\end{bmatrix}
\begin{bmatrix} u_1^j \ u_2^j \u_3^j \end{bmatrix}\
&=(v_1^{(j)},v_2^{(j)}, v_3^{(j)}) \cdot K_e \cdot \begin{bmatrix} u_1^j \ u_2^j \u_3^j \end{bmatrix}
\end{align}
$$

十二、面积坐标的梯度

对面积坐标 $L_i(x,y)$，通过对坐标变换矩阵求导，得梯度分量：
$$
\frac{\partial L_i}{\partial x} = \frac{y_j - y_k}{2S}, \quad \frac{\partial L_i}{\partial y} = \frac{x_k - x_j}{2S}
$$
$(i,j,k)$ 是三角形顶点的循环标号，$S$ 是三角形面积

性质：$\nabla L_i$ 是常向量（与位置无关 ）

十三、刚度矩阵的积分

单元刚度矩阵的积分项（梯度内积）可简化为：
$$
\int_{\Delta_j} \nabla L_k \cdot \nabla L_l , dxdy = \nabla L_k \cdot \nabla L_l \cdot S
$$
利用 $\nabla L_i$ 是常向量，积分转化为向量点乘与面积乘积

十四、单元载荷的积分

载荷项 $\int_{\Delta_j} f u_h , dxdy$ 分解为基函数的线性组合，积分后表示为：
$$
\begin{align}
\int_{\Delta_j} f u_h , dxdy &= \int_{\Delta_j} f (v_1^{(j)}L_1^{(j)}+v_2^{(j)}L_2^{(j)}+v_3^{(j)}L_3^{(j)}), dxdy \
&=(v_1^{(j)},v_2^{(j)},v_3^{(j)})\begin{bmatrix}\displaystyle \int_{\Delta_j} f\cdot L_1^{(j)},dxdy \ \displaystyle \int_{\Delta_j} f\cdot L_2^{(j)},dxdy \ \displaystyle \int_{\Delta_j} f\cdot L_3^{(j)},dxdy \end{bmatrix}
\end{align}
$$