3.FEM
两块无线大的PEC板位于YOZ平面，一块位于x=0电位为0V，另一块位于x=4电位为20V，使用一维有限元方法求解板间电位。
（1）求解电位的解析表达式
由泊松方程
两边积分2次得到
由边界条件
（2）Ritz变分法FEM
利用讲义（21）式
将区域离散化后，在每个子区域上的
代入（21）式得到
然后对
得到含有N-2个未知量的N-2个线性方程组，求解线性方程组即得到离散点上的解，然后由
编写Matlab程序如下：
clc; clear;
N=20;
xstart=0; xend=4; %x range
xx=linspace(xstart, xend, N);
syms x y F1 F2
%construct y
y=ones(1,N);
y=sym(y);
for i=1:N
y(i)=['y', num2str(i)];
end
y(1)=0; y(N)=20; %bound
f=x^2 + 1/2*x + 1/4;
tic
% Ritz method
for i=1:N-1
Fd(i) = 0.5*int( ( ( y(i+1)-y(i) ) / ( xx(i+1)-xx(i) ) )^2, x, xx(i), xx(i+1) )...
+int( f * ( y(i) * ( xx(i+1)-x ) / ( xx(i+1)-xx(i) ) + y(i+1) * ( x-xx(i) ) / ( xx(i+1)-xx(i) ) ) , x, xx(i), xx(i+1) );
end
F=sum(Fd);
for i=2:N-1
Fdiff(i)=collect(diff(F, y(i)));
for j=2:N-1
temp=coeffs(Fdiff(i), y(j)); % Extract the coefficient matrix A, Ax=b
if( size(temp) == 1 )
A(i-1,j-1) = 0;
else
A(i-1,j-1) = temp(2);
end
temp=coeffs(Fdiff(i)); % Extract the right matrix b
b(i-1)=temp(1);
end
end
A=double(A); b=double(b);
yy(2:N-1)=inv(A)*(-b');
t=toc
yy(1)=double( y(1) );
yy(N)=double( y(N) );
plot( xx, yy, 'b*'); hold on
plot( xx, yy,'b--');
%Analytic solution
ax=xstart:0.001:xend;
ay=1/12*ax.^4 + 1/12*ax.^3 + 1/8*ax.^2 - 13/6*ax;
plot(ax, ay, 'r')
title(['Ritz FEM 电位分布，N=', num2str(N)]);
xlabel('距离'); ylabel('电位');
legend('精确数值解','数值解', '解析解');
N=5,10,20,30时的电位分布分别如下图：
（3）Garlerkin法FEM
利用讲义（29）式，加权余量为：
将区域离散化后，在每个子区域上的
检验函数
将
令
得到含有N-2个未知量的N-2个线性方程组，求解线性方程组即得到离散点上的解，然后由
编写Matlab程序如下：
clc; clear;
N=20;
xstart=0; xend=4; %x range
xx=linspace(xstart, xend, N);
syms x F1 F2
%construct y
y=ones(1,N);
y=sym(y);
for i=1:N
y(i)=['y', num2str(i)];
end
y(1)=0; y(N)=20; %bound
f=x^2 + 1/2*x + 1/4;
tic
% Galerkin method
for i=2:N-1
F(i) = -int( ( ( y(i)-y(i-1) ) / ( xx(i)-xx(i-1) ) ) * 1 / ( xx(i)-xx(i-1) ), x, xx(i-1), xx(i) )...
-int( ( ( y(i+1)-y(i) ) / ( xx(i+1)-xx(i) ) ) * (-1) / ( xx(i+1)-xx(i) ), x, xx(i), xx(i+1) )...
-int( f * ( x-xx(i-1) ) / ( xx(i)-xx(i-1) ) , x, xx(i-1), xx(i) )...
-int( f * ( xx(i+1)-x ) / ( xx(i+1)-xx(i) ) , x, xx(i), xx(i+1) );
end
for i=2:N-1
Ff(i)=collect(F(i));
for j=2:N-1
temp=coeffs(Ff(i), y(j)); % Extract the coefficient matrix A, Ax=b
if( size(temp) == 1 )
A(i-1,j-1) = 0;
else
A(i-1,j-1) = temp(2);
end
temp=coeffs(Ff(i)); % Extract the right matrix b
b(i-1)=temp(1);
end
end
A=double(A); b=double(b);
yy(2:N-1)=inv(A)*(-b');
t=toc
yy(1)=double( y(1) );
yy(N)=double( y(N) );
plot( xx, yy, 'b*', xx, yy, 'b--'); hold on
%Analytic solution
ax=xstart:0.001:xend;
ay=1/12*ax.^4 + 1/12*ax.^3 + 1/8*ax.^2 - 13/6*ax;
plot(ax, ay, 'r')
title(['Galerkin FEM 电位分布，N=', num2str(N)]);
xlabel('距离'); ylabel('电位');
legend('精确数值解','数值解', '解析解');
N=5,10,20,30时的电位分布分别如下图：
（4）比较Ritz和Garlerkin法FEM的收敛性
两种方法均收敛。
收敛速度（使用Matlab的tic，toc计时积分、提取系数矩阵和解方程的时间）：
N=10时，Ritz用时0.266，Garlerkin用时0.297
N=20时，Ritz用时0.953，Garlerkin用时1.078
N=30时，Ritz用时2.047，Garlerkin用时2.219
所以，Ritz方法收敛速度要快于Garlerkin方法。