该博客详细介绍了如何利用Matlab解决二维对流方程，涉及A、B、C格式和蛙跳格式的差分方法。通过改进MeshGrid、ParasValidater、TryGetElement等函数，优化了代码效率。同时，博主分析了在不同时间步长下，四种格式的稳定性，并展示了结果动画。
1. 题目
2. 转述
原题目要求可以简化为：
对于二维对流方程：
∂u/∂t+∂u/∂x+∂u/∂y=0
u(x,y,0)={█(1,when-4≤x≤4,-4≤y≤4@0,other)┤
取计算范围为 -16≤x≤16,-16≤y≤16，Δx=Δy=1，Δt/Δx=Δt/Δy=0.5,1,2，t>0，使用 A，B，C格式与蛙跳格式计算90个时间步长，制作 u-x-y 图随时间的动画，依此分析四种格式对于二维对流方程的稳定性。
3. 分析
之前分析一维对流方程的代码
https://blog.csdn.net/PriceCheap/article/details/124172427
基本可以照抄的
只不过现在看还是有可以修改的地方
3.1 MeshGrid
之前求 X F 向量是封到一个函数里
旧代码 1
% 多个 dxdt 情况对比
for i = 1:1:size(dxdt,2)
Delta_t = Delta_x/dxdt(i);
[F,X,T] = PreDifferenceProcess(x_border,x_initialBorder,Delta_x,Delta_t,@InitialCondition);
旧代码 2
function [F,X,T] = PreDifferenceProcess(x_border,x_initialBorder,Delta_x,Delta_t,InitialCondition)
% 求差分格式之前需要的前处理
% t 的最大步数
% 由 x_initialBorder 与 x_border 之间的空隙决定
n_border = (x_border - x_initialBorder)/Delta_x + 1;
% x 的向量
X = -x_border : Delta_x : x_border;
% t 的向量
T = 0 : Delta_t : n_border * Delta_t;
% 流场物理量矩阵初始化
F = zeros(size(X,2),size(T,2));
% 初始条件
for i = 1:1:size(F,1)
F(i,1) = InitialCondition(X(i));
end
end
现在想想是没有必要的，直接 MeshGrid 就好了
这样直接用 X Y 矩阵也挺好
新代码
% 求 T 向量
t = differenceBorder(5):dt:differenceBorder(6);
% 多个 dxdt 情况对比
for i = 1:1:max(size(dtdx))
% 求步长
dx = dt/dtdx(i);
dy = dt/dtdy(i);
% 求 X 和 Y 向量
x = differenceBorder(1):dx:differenceBorder(2);
y = differenceBorder(3):dy:differenceBorder(4);
% 求 X 和 Y 矩阵
% X 矩阵每一行都是 X 向量
% Y 矩阵每一列都是 Y 向量
% 和 XOY 坐标系的直觉相配
[X,Y] = meshgrid(x,y);
3.2 ParasValidater
闲得无聊写了一个验证输入的函数……万一有人真的想用我的函数呢x
类似这样
% 验证 alpha, beta 是否合法
if alpha == 0
error("alpha 需为非零！");
end
if beta == 0
error("beta 需为非零！");
end
% 验证 differenceBorder 是否合法
if min(size(differenceBorder)) ~= 1 || length(size(differenceBorder)) ~= 2
error("differenceBorder 需为一维向量！");
end
if max(size(differenceBorder)) ~= 6
error("differenceBorder 需有六个元素！");
end
if differenceBorder(1) > differenceBorder(2)
error("differenceBorder 中 X 的下界需要小于 X 的上界！");
end
if differenceBorder(3) > differenceBorder(4)
error("differenceBorder 中 Y 的下界需要小于 Y 的上界！");
end
if differenceBorder(5) > differenceBorder(6)
error("differenceBorder 中 T 的下界需要小于 T 的上界！");
end
3.3 TryGetElement
差分格式要取的点可能超出边界
我之前的做法是生硬地分情况讨论
旧代码
if i<2
tmp2 = 0;
else
tmp2 = F(i-1,n);
end
if i>size(F,1)-1
tmp1 = 0;
else
tmp1 = F(i+1,n);
end
现在我统一使用一个函数来判断是否需要取 0
新代码 2
for x = 1:1:size(Zeta,1)
for y = 1:1:size(Zeta,2)
% 基点
base = TryGetElement(Zeta,x,y);
% 基点 x 正方向一步
right = TryGetElement(Zeta,x+1,y);
% 基点 x 负方向一步
left = TryGetElement(Zeta,x-1,y);
% 基点 y 正方向一步
up = TryGetElement(Zeta,x,y+1);
% 基点 y 负方向一步
down = TryGetElement(Zeta,x,y-1);
新代码 2
function [element] = TryGetElement(Array,row,col)
% 超出边界的可以取 0 而不损失精度
if row < 1 || row > size(Array,1)
element = 0;
return;
end
if col < 1 || col > size(Array,2)
element = 0;
return;
end
element = Array(row,col);
end
3.4 Zeta(x,y)
之前我是把待求函数的值和时间放在一起了，所以就显得很……乱吧
这下我把这两个分离开了，所以我可以写 Zeta(x,y) 而不是 Zeta(x,y,n+1)
旧代码 A 格式
function [F] = DifferenceFormat_A(F,a,Delta_x,Delta_t)
% A 格式
% (F(i,n+1)-F(i,n))/Delta_t + a*(F(i+1,n)-F(i-1,n))/(2*Delta_x) = 0
% => F(i,n+1) = F(i,n) -a/2*Delta_t/Delta_x*(F(i+1,n)-F(i-1,n))
for n = 1:1:size(F,2)-1
for i = 1:1:size(F,1)
% 超出边界的可以取 0 而不损失精度
if i<2
tmp2 = 0;
else
tmp2 = F(i-1,n);
end
if i>size(F,1)-1
tmp1 = 0;
else
tmp1 = F(i+1,n);
end
F(i,n+1) = F(i,n) - a/2*Delta_t/Delta_x*(tmp1-tmp2);
end
end
end
新代码 A 格式
function [Zeta_New] = DifferenceFormat_A_2D(Zeta,alpha,beta,dtdx,dtdy)
% 二维对流方程 A 格式
% (Zeta(xi,yi,n+1)-Zeta(xi,yi,n))/dt + alpha*(Zeta(xi+1,yi,n)-Zeta(xi,yi,n))/(2*dx) + beta*(Zeta(xi,yi+1,n)-Zeta(xi,yi,n))/(2*dy) = 0
% => Zeta(xi,yi,n+1) =
% Zeta(xi,yi,n) - alpha/2*dt/dx*(Zeta(xi+1,yi,n)-Zeta(xi-1,yi,n)) - beta/2*dt/dy*(Zeta(xi,yi+1,n)-Zeta(xi,yi-1,n))
Zeta_New = zeros(size(Zeta));
for x = 1:1:size(Zeta,1)
for y = 1:1:size(Zeta,2)
% 基点
base = TryGetElement(Zeta,x,y);
% 基点 x 正方向一步
right = TryGetElement(Zeta,x+1,y);
% 基点 x 负方向一步
left = TryGetElement(Zeta,x-1,y);
% 基点 y 正方向一步
up = TryGetElement(Zeta,x,y+1);
% 基点 y 负方向一步
down = TryGetElement(Zeta,x,y-1);
Zeta_New(x,y) = base - alpha/2*dtdx*(right-left) - beta/2*dtdy*(up-down);
end
end
end
3.5 dtdx dtdy
之前我在总的差分函数里面是传入一个 dtdx 的向量，然后我函数里面对这个向量 for
但是现在我有了打印 png 的需求，就感觉，什么地方都要用到 dtdx 的计数器 i 来区分情况就很烦
所以感觉把这个总的函数做成只处理 dtdx dtdy 为常数的就好了，这才是简洁的
在外面多次调用这个函数
旧代码 ODCE_DifferenceProcess
% 多个 dxdt 情况对比
for i = 1:1:size(dxdt,2)
Delta_t = Delta_x/dxdt(i);
[F] = DifferenceFormat(F,a,Delta_x,Delta_t);
新代码 ConvectionEqDifferenceMethod_2D
Zeta = DifferenceFormatFun(Zeta,alpha,beta,dtdx,dtdy);
旧代码 外部调用函数
ODCE_DifferenceProcess(@DifferenceFormat_A,'A',@InitialCondition,a,x_border,x_initialBorder,Delta_x,dxdt);
新代码 外部调用函数
for i = 1:1:max(size(dtdx))
ConvectionEqDifferenceMethod_2D(@DifferenceFormat_A_2D,'A', ...
alpha,beta,differenceBorder, ...
@InitializationFun,initializationBorder, ...
dt,dtdx(i),dtdy(i))
end
3.6 moive
创建动画的话，我本来还想使用 print 输出多个 png，然后再用 ps 拼成一个 gif 的
之后发现了 movie 可以直接播放动画，那就不用 gif 这么麻烦了
直接在每一次画图的时候先 drawnow 然后 getframe 最后播放就好了
旧代码 ODCE_DifferenceProcess
% 多个 dxdt 情况, 多个缩放情况对比
for i = 1:1:size(dxdt,2)
figure();
for j = 1:1:size(scale,2)
subplot(size(scale,2),1,j);
mesh(T_scaled,X_scaled,F_scaled,'FaceColor','flat');
新代码 ConvectionEqDifferenceMethod_2D
% 迭代时间步
for i = 1:1:max(size(t))
% ---画图---
% 平面着色
mesh(X,Y,Zeta,'FaceColor','flat');
% 捕获影片帧
drawnow;
Moive(i) = getframe;
3.7 InRange
题目中给的初始化函数可以简化为一个是否在范围内的函数
Matlab 没有提供这个函数我有点奇怪
或许是我查询的关键字错了才没有查到
function [IsInRange] = InRange(element,lower,upper)
% 在范围内返回 1
% 不在范围内返回 0
if element > lower && element < upper
IsInRange = 1;
return;
end
IsInRange = 0;
end
3.8 No Searching Loop
我在别的地方看到一个原则就是
推荐使用矩阵矢量化运算，而不是多重循环遍历矩阵
因为矩阵是 Matalb 的一个结构体，对矢量运算有优化
但是多重循环遍历矩阵每个元素那就一定是 o(nm) 的时间
我定义了一个 InRange 来判断一个函数是否是在一个给定的范围内
如果直接用的话，那么如果输入一个矩阵，这个矩阵会和一个数比较，比较不了
测试代码 1
InitializationFun = @(X,Y) InRange(X,initializationBorder(1),initializationBorder(2)).*InRange(Y,initializationBorder(3),initializationBorder(4));
我本来想用 bsxfun 来实现不用循环，因为它的功能是对两个矩阵的每一组对应元素都执行给定函数嘛
结果我发现即使如此，bsxfun 传给输入函数的变量依然是矩阵而不是元素
测试代码 2
InInitializationRange = @(x,y) InRange(x,initializationBorder(1),initializationBorder(
最低0.47元/天 解锁文章
4036
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
