二维波动方程的简单数值解（Matlab）
贡献者： addis
- 本文处于草稿阶段。
图 1：波的干涉，动画见这里。
\begin{equation}
\boldsymbol{\nabla}^2 u - \frac{1}{v_0^2} \frac{\partial^{2}{u}}{\partial{t}^{2}} = 0~,
\end{equation}
\begin{equation}
\frac{\partial^{2}{u}}{\partial{t}^{2}} = v_0^2 \boldsymbol{\nabla}^2 u~.
\end{equation}
使用差分法 求解常微分方程，时间差分得
\begin{equation}
u(t_{i+1}) = v_0^2\Delta t^2 \boldsymbol{\nabla}^2 u(t_i) + 2u(t_{i}) - u(t_{i-1})~.
\end{equation}
\begin{equation}
\boldsymbol{\nabla}^2 u(x_i,y_i) = \frac{u(x_{i+1}, y_i) + u(x_{i-1}, y_i) - u(x_{i}, y_i)}{\Delta x^2}
+ \frac{u(x_i, y_{i+1}) + u(x_i, y_{i-1}) - u(x_i, y_{i})}{\Delta y^2}~.
\end{equation}
边界条件：反射，行波边界条件（类比式 5 ）：
对平面波 $u = A \sin\left( \boldsymbol{\mathbf{k}} \boldsymbol\cdot \boldsymbol{\mathbf{r}} - \omega t\right) $，有
\begin{equation}
v_0 \boldsymbol\nabla u = - \hat{\boldsymbol{\mathbf{k}}} \frac{\partial z}{\partial t} ~,
\end{equation}
\begin{equation}
\frac{\partial z}{\partial t} = -s v_0 \left\lvert \boldsymbol\nabla u \right\rvert ~.
\end{equation}
以下 Matlab 代码中使用了
代码 1：wav2D.m
% 二维波动方程的简单数值解
% 使用简单的差分法
function wav2D
% === 参数设置 ====
bc = 'o'; % 边界条件: [d] Dirichlet, [n] Neumann, [o] Open
v0 = 1;
Nplot = 40;
tmin = 0; tmax = 20; Nt = 4001;
xmin = -10; xmax = 10; Nx = 800;
ymin = -10; ymax = 10; Ny = 800;
crange = [-0.4,0.4];
% ================
close all;
t = linspace(tmin, tmax, Nt); dt = (tmax-tmin)/(Nt-1);
x = linspace(xmin, xmax, Nx); dx = (xmax-xmin)/(Nx-1);
y = linspace(ymin, ymax, Ny); dy = (ymax-ymin)/(Ny-1);
% [X, Y] = ndgrid(x, y);
Z = zeros(Nx, Ny, 3);
figure; set(gcf, 'Unit', 'Normalized', 'Position', [0.1,0.1,0.4,0.4]);
for it = 1:Nt
Z(400,400,2) = 2*sin(5*t(it));
d2Z = laplacian(Z(:,:,2), dx, dy);
if bc == 'd'
Z(:,:,3) = v0^2*dt^2*d2Z + 2*Z(:,:,2) - Z(:,:,1);
elseif bc == 'o'
Z(2:end-1,2:end-1,3) = (v0*dt)^2*d2Z(2:end-1, 2:end-1) +...
2*Z(2:end-1,2:end-1,2) - Z(2:end-1,2:end-1,1);
% ymax
grad_x = gradient(Z(:,end,2))/dx;
grad_y = (Z(:,end,2)-Z(:,end-1,2))/dy;
Z(:,end,3) = Z(:,end,2) ...
- (v0*dt)*sign(grad_y).*sqrt(grad_x.^2 + grad_y.^2);
% ymin
grad_x = gradient(Z(:,1,2))/dx;
grad_y = (Z(:,2,2)-Z(:,1,2))/dy;
Z(:,1,3) = Z(:,1,2) ...
+ (v0*dt)*sign(grad_y).*sqrt(grad_x.^2 + grad_y.^2);
% xmax
grad_x = (Z(end,:,2)-Z(end-1,:,2))/dx;
grad_y = gradient(Z(end,:,2))/dy;
Z(end,:,3) = Z(end,:,2) ...
- (v0*dt)*sign(grad_x).*sqrt(grad_x.^2 + grad_y.^2);
% xmin
grad_x = (Z(2,:,2)-Z(1,:,2))/dx;
grad_y = gradient(Z(1,:,2))/dy;
Z(1,:,3) = Z(1,:,2) ...
+ (v0*dt)*sign(grad_x).*sqrt(grad_x.^2 + grad_y.^2);
end
Z(:,:,1) = Z(:,:,2); Z(:,:,2) = Z(:,:,3);
% 画图
if mod(it, Nplot) == 0
clf; surfCart(x, y, Z(:,:,3)); caxis(crange);
colormap jet; shading interp; axis equal;
xlabel x; ylabel y; view(78, 32);
axis([xmin,xmax,ymin,ymax,-2,2]);
title(['t = ' num2str(t(it), '%.2f')]);
saveas(gcf, [num2str(it/Nplot) '.png']);
end
end
end
% laplacian of Z(ix, iy)
% finite difference
% boundary condition (not included in Z) is 0
function d2Z = laplacian(Z, dx, dy)
N1 = size(Z, 1); N2 = size(Z, 2);
z1 = zeros(N1, 1); z2 = zeros(1, N2);
d2Z_1 = ([Z(2:end, :); z2] + [z2; Z(1:end-1, :)] - 2*Z)/dx^2;
d2Z_2 = ([Z(:, 2:end) z1] + [z1 Z(:, 1:end-1)] - 2*Z)/dy^2;
d2Z = d2Z_1 + d2Z_2;
end
未完成：这个代码应该放到别处并引用
代码 2：surfCart.m
% my version of surf()
% X, Y, Z are 2D matrices
% the color of each facet is determined by the value
% of the upper-left grid point
% there will be 1 facet for each element of Z,
% and one more grid-point at
% the end of each dimension
function h = surfCart(x,y,Z,varargin)
[X,Y] = ndgrid([x,2*x(end)-x(end-1)],[y,2*y(end)-y(end-1)]);
Z = [Z, Z(:,end); Z(end,:), Z(end,end)];
h = surf(X,Y,Z,varargin{:});
shading flat;
view(90,90);
axis([min(X(:)), max(X(:)), min(Y(:)), max(Y(:)), ...
[-1, 1]*max(1, max(abs(Z(:))))]);
xlabel x; ylabel y; zlabel z;
colorbar;
set(datacursormode(gcf), 'UpdateFcn', @datatip);
cameratoolbar;
end