在 Matlab 上的案例 P2 中，我致力于使用 FEM 解决 PDE。我尝试使用二次拉格朗日形状函数正确组装负载矢量

[rspts,qwgts] = Gausspoints(precision); % quadrature rule
np      = size(p,2);                   % number of nodes
nt      = size(t,2);                   % number of elements

b = sparse(np,1);

for i=1:nt % loop over elements

nodes = t(1:6,i); % node numbers
x = p(1,nodes);   % node x-coordinates
y = p(2,nodes);   % node y-coordinates

bK=zeros(6,1); % elements load vector

for q=1:length(qwgts)                                  % quadrature loop   
r=rspts(q,1);                                         % quadrature r-coordinate
s=rspts(q,2);                                         % quadrature s-coordinate
[S,dSdx,dSdy,detJ]=Isopmap(x,y,r,s,@P2shapes);        % map

wxarea=qwgts(q)*detJ/2;                               % weight times area

bK=bK+S*f(mean(r),mean(s))*wxarea; % elements load vector
end
b(nodes)       = b(nodes) +bK;
end
   

我的问题... 当我采用源项 f=0 的 PDE 时，数值解和精确解的图给出了相同的结果。但是，如果我输入例如 f(x,y)=x+y，我会以两个解决方案之间的情节差异结束。

简要说明 我找到了有用的信息：计算可变系数的刚度矩阵

解决方案： 我没有正确理解将（x，y）域映射到（r，s）域的逻辑。积分可以这样写：（正好在元素负载向量的行中）

 ...

x_physical= dot(x,S);
y_physical= dot(y,S);
bK=bK+S*f(x_physical,y_physical)'*wxarea;
...