萌新学习有限元方法解偏微分方程编程实例（一）

本人新手入门，从最简单最基础的例子学习，不足之处请大家多多批评指正。

以何晓明老师的课件例题为例：例题一

运用MATLAB编写驱动文件run.m和函数文件unit.m

%% 初始化
clc,clear
N=4;
A=zeros(N+1,N+1);
B=zeros(N+1,1);
xx=[0:1:N];
h=1.0/N;
xxx=xx.*h;

%真实解
uu = @(x) x.*cos(x);

%% 逐个循环，生成单元刚度矩阵，并且组装总刚度矩阵
for j = 1:N
   x1=xxx(j);x2=xxx(j+1); 
   [AA,BB] = unit(N,h,j,x1,x2);
   A(j:j+1,j:j+1)=AA+A(j:j+1,j:j+1);
   B(j:j+1,  1  )=BB+B(j:j+1,  1  );
end

%% 边界条件
   A(1,:)=0;
   A(1,1)=1;
   A(N+1,:)=0;
   A(N+1,N+1)=1; 
   B(1)=uu(0);
   B(N+1)=uu(1);
  
%% 求解,作图,误差
u = inv(A) * B;
uu(xxx);
plot(xxx,u,'b--',xxx,uu(xxx),'r-')
legend('数值解(蓝色虚线)','真实解(红色实线)','Location', 'best');
error = max(u-uu(xxx)');
 

function [AA,BB] = unit(N,h,j,x1,x2)

% unit.m生成单个单元刚度矩阵
% 请修改为自己的c(x) 
% AA是单元刚度矩阵
% FF是单元载荷矩阵

c =  @(x) exp(x)*(1/h)*(1/h);
FB1 = @(x) -exp(x).*(cos(x)-2.*sin(x)-x.*cos(x)-x.*sin(x)).*(x2-x)*(1/h);
FB2 = @(x) -exp(x).*(cos(x)-2.*sin(x)-x.*cos(x)-x.*sin(x)).*(x-x1)*(1/h); 
 
AA=zeros(2,2);
BB=zeros(2,1);
AA(1,1)=  integral(c,x1,x2);
AA(1,2)= -integral(c,x1,x2);
AA(2,1)= -integral(c,x1,x2);
AA(2,2)=  integral(c,x1,x2);

BB(1,1)=  integral(FB1,x1,x2);
BB(2,1)=  integral(FB2,x1,x2);

end

古典解和数值解的拟合图像如下（N=4）：

 可以验证误差（error）满足下列表格的数值：