有限差分方法详细思路和步骤+Matlab高效编程求解技巧

有限差分法在（偏）微分方程数值解方面可谓是举足轻重。我们在学习导数的概念时知道，连续函数的导数其实是通过相邻离散点作差分然后取极限得到的，即：

解微分方程时，将所有的导数项用离散点处的差分来近似，这就是有限差分法的主要思想。选择不同的离散点来将    处的导数    作近似，就会得到不同的有限差分格式。比如，上面的式子我们是用    点与    点做差分，这种叫前向差分(forward difference)；你也可以用    点与    点做差来近似    ：

这叫后向差分(backward difference)；还有许多其他的格式，不同的格式都是在追求同一个目标：精度更好求解更稳定，也就是相同离散点的情况下能更准确地对导数进行近似。

一阶导数简单，那二阶导数甚至高阶导数呢？其实也很简单，因为，根据定义，你一定知道二阶导数是导数的导数：

(如果    那就是等号）

紧接着我们将分子中的一阶导数再用一次前向有限差分：

因此立刻可以写出二阶导数差分公式：

也就是说，二阶导数可由    点及其左右相邻点一共3个离散点处的函数值来近似（你用泰勒公式展开，将步长h看作小量，肯定能得到    ，不信试试哈）

同样地，不同的近似算法对应不同的差分格式，上述对应最常用的中心差分格式，还有精度更高的五点差分格式等，在此不表。怎样的差分格式不重要，重要的是你要能弄明白到底用到哪些离散点来计算高阶导数，这样方便你在编程中体现。

用实例手把手教学

下面我们用一个简单的波动方程来手把手教会读者：

• 如何利用有限差分方法进行离散化
• 如何写简洁高效的MATLAB程序
• 如何进行数值解的可视化

琴弦振动问题

依然用PINN神经网络求解PDE一文中的方程：

边界条件

初始条件

这是一个具有明确物理意义的问题：一条琴弦，两端固定在0和1的位置，拉动琴弦使得它初始具有构型    ：

那么接下来琴弦肯定要振动对吧，那它将如何振动呢？这就是上面的方程干的事。

连续方程的有限差分离散

时间、空间离散化，初始边界条件如何处理、方程的离散格式如何给出内部未知场函数，详细可见👇

结果的可视化

基于上面的离散求解方法，我们得到琴弦振动动态图：

动图的生成，只要利用一个for循环，搭配上drawnow指令即可实现，推荐读者查看往期推文如何在matlab中生成动图视频

简介高效的MATLAB代码

下面给出完整求解代码：

c=2;
deltax=0.01;
deltat=deltax/5;
Nx=1/deltax+1; % total points along x
x=0:deltax:1;
T=1.5; % total time
Nt=round(T/deltat+1); % total points along t

U=zeros(Nt,Nx);
U(1,:)=3*x.*(x<=1/3)+3/2*(1-x).*(x>=1/3); % the first IC
U(2,:)=U(1,:); % the second IC  
U(:,1)=0;
U(:,end)=0; % BCs

r=c^2*(deltat/deltax)^2;
for n=2:Nt
    U(n+1,2:end-1)=-U(n-1,2:end-1)+r*U(n,1:end-2)+(2-2*r)*U(n,2:end-1)+r*U(n,3:end); % 离散格式
end

读者需要注意或者学习的是最后for循环里面的一句长指令，这句指令将看似比较麻烦的离散方程

高效地用代码表述，相信我，这是非常有用的Matlab高效编程技巧！

这是生成动图的代码：

figure
for i=1:Nt
    plot(x,U(i,:),'LineWidth',1.4)
    xlabel('x')
    ylabel('u')
    ylim([-1.5,1.5])
    set(gca,'PlotBoxAspectRatio',[1 0.4 1],'FontSize',15)
    title(['t=',num2str((i-1)*deltat)])
    drawnow
end

以后你遇到要生成可视化动图的问题，都可以借鉴这段代码。