MATLAB绘制球谐函数报错：网格数组需为NDGRID结构

问题背景 #

我在3D声学与流体流动场景中，想用球谐函数表示球体的扰动。目前已实现3D球体的xy平面2D扰动，现在要扩展到包含z方向的扰动。

现有可运行的2D扰动代码 #

clearvars; clc; close all;
Nx = 128;
Ny = 128;
Nz = 128; 
Lx =128; 
Ly = 128; 
Lz = 128; 
xi  = (0:Nx-1)/Nx*2*pi;
xi_x =  2*pi/Lx;
x =  xi/xi_x;
yi  = (0:Ny-1)/Ny*2*pi;
yi_y =  2*pi/Ly;
y = yi/yi_y;
zi  = (0:Nz-1)/Nz*2*pi;
zi_z =  2*pi/Lz;
z = zi/zi_z;
[X,Y,Z] = meshgrid(x,y,z);
A = 2*pi / Lx;
B = 2*pi / Ly;
C = 2*pi / Lz;
x0 = 64; 
y0 = 64;
z0 = 64; 
rx0 = 20; 
ry0 =  20; 
rz0 =  20;
p = 3;
b = 0.1; % pert amplitude 
c = 12; 
d = 1;
a = 4;
theta = atan2(Y -y0,  X-x0) - (pi/c);
p0 = ((X-x0) .* (X-x0)) /(rx0 * rx0) + ((Y-y0) .* (Y-y0))/(ry0 * ry0) + ((Z-z0) .* (Z-z0))/(rz0 * rz0);
Test =d + a * exp((-1. * p0 .* (1 - b .* cos(c * theta))).^p) ;
figure
isosurface(X,Y,Z,Test);
shading flat;
grid on;

尝试的球谐函数扩展代码 #

clearvars; clc; close all;
%in spherical coord
%calculate r
Nx = 128;  
Ny = 128;  
Nz = 128;
Lx =128; 
Ly = 128; 
Lz = 128; 
xi  = (0:Nx-1)/Nx*2*pi;
xi_x =  2*pi/Lx;
x =  xi/xi_x;
yi  = (0:Ny-1)/Ny*2*pi;
yi_y =  2*pi/Ly;
y = yi/yi_y;
zi  = (0:Nz-1)/Nz*2*pi;
zi_z =  2*pi/Lz;
z = zi/zi_z;
r = sqrt(x.^2 + y.^2 + z.^2); 
% Create the grid
delta = pi/127; 
%Taking for instance l=1, m=-1 you can generate this harmonic on a (azimuth, elevation) grid like this:
azimuths = 0 : delta : pi; 
elevations = 0 : 2*delta : 2*pi; 
[R, A, E] = ndgrid(r, azimuths, elevations); %A is phi and E is theta
H = 0.25 * sqrt(3/(2*pi)) .* exp(-1j*A) .* sin(E) .* cos(E);
%transform the grid back to cartesian grid like this:
%can also add some radial distortion to make things look nicer:
%the radial part depends on your domain
X = r .* cos(A) .* sin(E); 
Y = r .* sin(A) .* sin(E); 
Z = r .* cos(E); 
%parameters
x0 = 64;  
y0 = 64; 
z0 = 64;  
rx0 = 20; 
ry0 =  20; 
rz0 =  20;
p = 3;
b = 0.1; % pert amplitude 
%c = 12; 
d = 1;
a = 4;
p0 = ((X-x0) .* (X-x0)) /(rx0 * rx0) + ((Y-y0) .* (Y-y0))/(ry0 * ry0) + ((Z-z0) .* (Z-z0))/(rz0 * rz0);
Test1 =d + a * exp((-1. * p0 .*H).^p) ;
figure
isosurface(X,Y,Z,real(Test1)); %ERROR

运行错误信息 #

Error using griddedInterpolant
Grid arrays must have NDGRID structure.

问题分析与解决 #

这个错误既不是球谐函数的设置问题，也不是Test1的函数形式问题，核心原因是网格结构不匹配：

• isosurface要求输入的X/Y/Z是meshgrid生成的三维笛卡尔网格（维度对应[Nx, Ny, Nz]），但你通过ndgrid(r, azimuths, elevations)转换后的X/Y/Z维度为[Nx, N_azimuth, N_elevation]，结构完全不符合要求。
• 额外问题：你用全局坐标计算球谐函数，而非以球心为原点的局部球面坐标，会导致扰动位置偏离目标球体。

修正后的代码示例 #

clearvars; clc; close all;
Nx = 128; Ny = 128; Nz = 128; 
Lx =128; Ly = 128; Lz = 128; 

% 生成标准笛卡尔网格（与原代码一致）
xi = (0:Nx-1)/Nx*2*pi; xi_x = 2*pi/Lx; x = xi/xi_x;
yi = (0:Ny-1)/Ny*2*pi; yi_y = 2*pi/Ly; y = yi/yi_y;
zi = (0:Nz-1)/Nz*2*pi; zi_z = 2*pi/Lz; z = zi/zi_z;
[X,Y,Z] = meshgrid(x,y,z);

% 球心参数
x0 = 64; y0 = 64; z0 = 64; 
rx0 = 20; ry0 = 20; rz0 = 20;

% 转换为以球心为原点的局部球面坐标
dx = X - x0; dy = Y - y0; dz = Z - z0;
r_local = sqrt(dx.^2 + dy.^2 + dz.^2);
theta = atan2(dy, dx); % 方位角φ
phi = acos(dz ./ r_local); % 极角θ（适配球谐函数定义）

% 计算球谐函数（取l=1,m=-1的实部作为z方向扰动）
harmonic = sqrt(3/(8*pi)) .* sin(phi) .* cos(theta);

% 叠加球谐扰动到原球体函数
p = 3; b = 0.1; d = 1; a = 4;
p0 = (dx.^2)/(rx0^2) + (dy.^2)/(ry0^2) + (dz.^2)/(rz0^2);
Test1 = d + a * exp( (-p0 .* (1 + b * harmonic)).^p );

% 绘制等值面
figure
isosurface(X,Y,Z,Test1);
shading flat; grid on; axis equal;

关键修正点 #

1. 保持笛卡尔网格结构与原代码一致，满足isosurface的输入要求
2. 基于球心局部坐标计算球面角度，确保扰动围绕目标球体
3. 取球谐函数的实部作为实数扰动，避免复数运算导致的异常

内容的提问来源于stack exchange，提问作者Jamie