[Matlab科学计算] 有限元法求二阶常微分方程

我们常常用微分方程来描述现实世界中的一些物理现象。由于微分方程的复杂性，只有在很简单的情况下才能得到微分方程的解析解。由于计算机的发展，采用数值方法求解微分方程的数值近似解得到广泛应用。微分方程的数值解法主要包括两大类：有限差分法和有限单元法。这里主要介绍有限单元法。

However，对于一个只学过微积分和矩阵论的工科生来说，要了解有限元法的数学原理还是有些困难，所以这里重点是介绍有限元法的思想和具体计算方法，深层的数学原理并不涉及，主要是本人也不懂~~~

1.微分方程基本概念

微分方程主要分为常微分方程和偏微分方程。
常微分方程：未知函数是单一自变量的函数，按照方程中导数的阶数，分为一阶、二阶、n阶常微分方程。
偏微分方程：未知函数是多个自变量的函数，至少是两个自变量，同样可以分为一阶、二阶、n阶偏微分微分方程。

二阶偏微分方程研究较多，又可以分为椭圆型方程、双曲型方程、抛物型方程等。这里不展开，具体可以查看偏微分方程数值解法相关书籍。

2.二阶常微分方程的有限元法求解举例

由于下面涉及大量公式，所以采用截图的方式进行介绍

3.具体Matlab实现如下：

%% main program
%一次有限元求一维常微分方程，基函数为分片线性函数
n = 10;
err = zeros(8, 1);% Linear interpolation
index = zeros(8, 1);
for i=1:8index(i) = 2^(i-1) * n; %10,20,40,80,160,320,640,1280
end
for i = 1:8N = index(i); %网格数量h = 1 / N; %网格大小x = 0:h:1; %节点向量xx = 0:1/1000:1;exact_solu = sin(pi.*xx); %精确解K = zeros(N-1); %初始化系数矩阵（刚度矩阵）K = K + diag(ones(N-1,1).*(2/h+2/3*h), 0) + diag((h/6-1/h).*ones(N-2,1), 1)...+ diag((h/6-1/h).*ones(N-2,1), -1);F = (1+pi^2)/(h*pi^2) * (2.*sin(pi.*x(2:N))-sin(pi.*x(1:N-1))-sin(pi.*x(3:N+1)))';K = sparse(K);coeff = K\F; %Ax=b,matlab求法x=A\bsolu = [0; coeff; 0]; %u(x)节点解%     figure()subplot(3,3,i)plot(x', solu, 'o--',xx', exact_solu, 'g--');legend('numerical solution', 'exact solution')xlabel('x')ylabel('U');xlim([-0.05,1.05])ylim([-0.05,1.05]);title(['numerical solution and exact solution N=',num2str(index(i))]);% calculate the errordu = solu(2:N+1) - solu(1:N);y1 = sin(2.*pi.*x);dy1 = y1(2:N+1) - y1(1:N);y2 = sin(pi.*x);dy2 = y2(2:N+1) - y2(1:N);error = pi^2*h/2*ones(1,N) + du'.^2/h + pi/4*dy1 - 2/h*dy2.*du';err(i) = sqrt(sum(error));
end
% figure()
subplot(3,3,9)
plot(log2(err),'o--');
title('误差阶')
xlabel('N')
ylabel('$log_{2}Error$', 'Interpreter','LaTex')
xlim([0.8,8.2])axis on;
set(gca,'xtick',1:1:8);
set(gca,'xticklabel',{5,10,20,40,80,160,320,640});[order,~] = polyfit(1:1:(length(err)), log2(err)', 1);
disp(order)

最后计算结果如下图所示：

参考文献：

代码 https://blog.csdn.net/waitingwinter/article/details/106164350