有限差分法求解含有高阶导数的微分方程：核心只要1行Matlab代码！
24-01-16 09:25发表于安徽
趁热打铁，本文向读者介绍含有高阶导数的微分方程如何利用有限差分法离散化，并提供简洁高效的matlab程序(十几行程序只有1行是核心)以及动态可视化图的代码👇
含有高阶导数的典型方程
本文考虑如下四阶微分方程
读者容易知道，此方程对应的物理问题是如下的两端固支的细长梁结构在均布载荷 作用下的横振动问题：
（ 在边界处的导数为0，意思就表明边界处转角为0，所以固支）
高阶导数的有限差分格式
在数值计算人的眼里，空间x变量不再是连续变量，而是一些均匀离散点，如下
根据有限差分法详细思路和步骤 ， 格点处的二阶空间导数写成这样（以 为中心的三点差分）：
这个没问题，但是上面方程中的 关于 的四阶导数该如何离散呢？你可以在心里把 记作 ，那么 的四阶导其实就是 的二阶导，而 在 处的二阶导为👇
因此，得到 的四阶导在 处为👇
再把 、 和 处二阶导的差分格式代入，即得到👇
也就是说，四阶导要用到 点及其左右一共五个离散点处的函数值，因此也叫五点差分公式。任意一个离散点 ，你只要利用它以及它左边两个点右边两个点共5个点即可表达出 处的四阶导数值。
用实例手把手实现
下面我们用上述梁横振动方程
来教读者：
如何利用有限差分方法进行离散化 如何写出简洁高效的MATLAB程序 如何进行数值解的动态可视化（振动动图）
连续方程的有限差分离散
首先依然要将时间和空间离散化👇
欲求解的场函数 即为一个 行、 列的数组。其中：
每一行代表一个时刻。因此，初始 时刻 意味着第一行为0；初始时刻 对t的导数也为0,意味着第二行也为0（为什么呢？因为初始时刻对t的导数可以写成 的差分格式）. 每一列代表一个x离散点.因此，边界 处 为0意味着第一列为0，边界 处 对x的导数为0，依然意味着第二列为0（为什么呢？因为x=0处y关于x的导数可以写成 的差分格式).同理， 处的边界条件可以知道倒数第一列和倒数第二列都为0.
下面再看中间离散点如何处理👇
也就是说，如果前面第 行和第 行已经知晓（比如初始条件给出了前两行值），那么第 的场函数值 就可以根据差分离散方程求出。要注意的是， 位置的空间四阶导数要用到左边两个点和右边两个点，因此这里的 肯定只能从第3列取，一直到倒数第3列为止。那前两列和最后两列咋办？看上面说的，边界条件已经把前两列和最后两列的值给出了😄✔
结果的可视化
基于上面的离散求解方法，我们得到梁的挠度 图：
以及梁的横向受迫振动动态图：
动图的生成，只要利用一个for循环，搭配上drawnow指令即可实现，推荐读者查看往期推文如何在matlab中生成动图视频或者高级跨年动画的制作
简洁高效的MATLAB代码
下面给出完整求解代码：
clc;close all
rho=7954;A=0.0336;I=0.00004032;E=5460000;
x=linspace(0,0.12,60);dx=x(2)-x(1);
t=linspace(0,0.1,80001);dt=t(2)-t(1); % 空间x和时间t离散点
nx=length(x);nt=length(t);
U=zeros(nt,nx);
U(1,:)=0; % 如果梁初始静止状态
U(2,:)=U(1,:); % 初始导数为0
U(:,1)=0; % 边界条件,x=0时y为0
U(:,end)=0; % 边界条件，x=0.12时，y为0
f=@(x,t)10*sin(60*t); % 右边加载函数
for n=2:nt-1
U(n+1,3:end-2)=1/(rho*A/dt^2)*(-rho*A/dt^2*(U(n-1,3:end-2)-2*U(n,3:end-2))+f(x,t(n))-E*I/dx^4*(U(n,1:end-4)-4*U(n,2:end-3)+6*U(n,3:end-2)-4*U(n,4:end-1)+U(n,5:end))); % 差分迭代下一行
U(n+1,2)=(3*U(n+1,1)+U(n+1,3))/4; % 由x=0导数边界条件，计算n+1行的第二个值
U(n+1,end-1)=(3*U(n+1,end)+U(n+1,end-2))/4; % 由x=0.12导数边界条件，计算n+1行的倒数第二个值
end
区区16行代码即搞定了这个问题！但是核心代码其实只有1句！那就是最后for循环里面的一句长指令。
这句指令将看似比较麻烦的离散差分方程仅用一行代码高效表达。具体来说，离散方程中我们要求解 ，如果你对 从3到 进行循环，那么效率就很low了。这里，我们用matlab的矩阵操作技巧，一次性地计算完成 到 的所有列元素，右边项中， 可以用U(n,4:end-1)一次性取出第n行的第4列到倒数第2列所有元素， 可以用U(n,5:end)一次性取出第n行的第5列到最后一列的所有元素,其他类似。笔者认为这是非常有用的Matlab高效编程技巧！任何一个想要提高matlab编程水平的读者都应掌握这种”向量化“编程的技巧。（读者不妨试一下如何用一句指令计算 在0到1上的二阶导数值）
最后贴上生成曲面图和动态图的代码：
% 结果可视化
figure
[X,T]=meshgrid(x,t);
surf(X,T,U)
shading interp
xlabel('x')
ylabel('t')
zlabel('y(x,t)')
% 动态图
figure
set(gcf,'Position',[400400700250],'Color',[111])
for ti=0:0.00025:t(end)
n=round(ti/dt)+1;
plot(x,U(n,:),'r','LineWidth',2)
xlabel('x')
title(['梁挠度随时间演化图，t=',num2str(ti)])
ylim([min(min(U)),max(max(U))])
drawnow
end
