MATLAB之物理场可视化

物理场指一个物理量的值在时间和空间中的分布，其定义是一个向量到另一个向量或数的映射。物理场分标量场和矢量场，对于标量场，在二维情况下可以使用等值线表示，在三维情况下则可以使用MATLAB绘图加强篇中提到的四维数据绘图方法，即使用颜色表示，当然也可以使用等值曲面表示；对于矢量场，我们则可以使用箭头、流线等来表示。

注：在矢量场中如果有需要也可根据场强大小绘制出等值线或等值面。

二维标量场

二维标量场使用等值线表示，函数contour可以绘制平面等值线。

z=peaks;		% peaks为MATLAB自带的演示函数
contour(peaks)		% 使用contourf函数绘制等值线

如果您觉得这样不够直观，那还可以使用contour3函数绘制立体等值线。(图片清晰度被降低了，读者可以自行去MATLAB跑一下代码，只需将contour替换为contour3即可)

又或者使用contourf函数在等值线内填充颜色。

三维标量场

四维数据作图

利用四维数据作图的技巧在MATLAB绘图加强篇已提到，此处直接给出代码，不再赘述。

clear all;clc;close all
[x,y,z]=meshgrid(-2:.2:2);	% .2=0.2,，meshgrid生成三维空间网格
v=x.*exp(-x.^2-y.^2-z.^2);	% 生成物理场
slice(v,[5 15],15,10)			% v所对应物理场的切片，分别在x=5,x=15,y=15,z=10处切片
axis([0 21 0 21 0 21]);		% 设置xyz坐标范围
hold on
colorbar('horiz')	% 在竖直方向显示温度表
colorbar('vert')		% 在水平方向显示温度表
view([-25 65])		% 通过方位角和俯仰角，设置观察角度

等值面

除了使用第四维数据作图来清晰的表示物理场之外，我们也可以仿照二维标量场，使用等值面来描绘物理场。

clc;clear all;close all
[x y z v]=flow;     % 数据来自MATLAB自带的演示函数
p=patch(isosurface(x,y,z,v,-3));    % 光照函数
isonormals(x,y,z,v,p)   % 数据归一化
set(p,'FaceColor','red','EdgeColor','none');    % 设置表面与边界颜色
daspect([1 1 1])    % 设置坐标轴比例
view(3)     % 视角
axis tight; grid on
camlight;   lighting phong  % 光照
alpha(.5)   % 设置表面透明度

二维矢量场

不同于二维标量场，矢量场在某点处既有大小又有方向，自然不能继续使用等值线表示，而是要使用 "带有箭头的等值线" ——流线以及大量的小箭头 来表示。

小箭头和流线的区别在于：小箭头只会根据某点的场强大小和方向两个属性绘图，小箭头仅表示该点的属性，与周围的场变化无关；而流线则不止受一个点的场强大小及方向影响，而流线密度同样反应了场的大小，因此流线可以更直观地反映出场的情况。

二维矢量场的箭头表示

二维矢量场可以使用quiver函数构成用来表示矢量场大小和方向的箭头，quiver(x,y,u,v)中x,y表示自变量，u,v表示矢量场的两个分量。

我们以电偶极子为例，假设有电量为 $4\pi \varepsilon_0$ 的正负电荷分别放置在 $(0.2, 0)$ 和 $(-0.2, 0)$ 处，形成一个电偶极子，则空间电势为 $V=\frac{1}{\sqrt{(x-0.2)^{2}+\sqrt{y-0.2)}2}}-\frac{1}{\sqrt{(x+0.2)^{2}+\sqrt{y-0.2)}2}}$，有了电势以后我们就可以求出其等势线和电场强度，进而将矢量场画出。

clc;clear all;close all
h=0.041; x=-0.3:h:0.3; y=-0.3:h:0.3;    % 此处h设置的是绘图精度，单独使用变量表示是为了方便后续调试代码
[X Y]=meshgrid(x,y);
Z=1./sqrt((X-0.2).^2+Y.^2)-1./sqrt((X+0.2).^2+Y.^2);	
[PX PY]=gradient(-Z,h);	% 利用梯度生成矢量场，电势的梯度为电场
contour(x,y,Z,[-12,-8,-5,-3,-2,-1,-.5,-0.1,0.1,0.5,1,2,3,5,8,12],'b')		% 绘制电势标量场的等值线——等势线
hold on
quiver(X,Y,PX,PY,'k')

二维矢量场的流线表示

streamline函数则可以构成流线，其所需的三组参数为空间坐标、矢量场的分量、流线起点。

clc;clear all;close all
load wind   % 数据为MATLAB内部的大气某处风速的数据
zmax=max(z(:));   zmin=min(z(:));
streamslice(x,y,z,u,v,w,[],[],(zmax-zmin)/2)
axis([70.1879, 134.3, 17.4999, 60]) 

三维矢量场

MATLAB提供了多种三维矢量场的可视化方法：可以用箭头quiver3或三维椎体coneplot表示各点场的大小；流线streamline，流线的切线表示该点场的大小，流线的密度表示场的大小；流管streamtube，流管的粗细反应矢量场的散度；流带streamribbon表示矢量场的旋转情况。

三维箭头 quiver3

clc;clear all;close all
t=0:0.5:8;  
x=sin(t);   y=cos(t);   z=t;
plot3(x,y,z,'linewidth',2)
hold on
u=gradient(x);  v=gradient(y);  w=gradient(z);
quiver3(x,y,z,u,v,w,0,'linewidth',2)
view(-60,60)

椎体 coneplot

clc;clear all;close all
load wind

% 确定要用于放置切片平面和指定圆锥体绘图所在位置的数据范围。
xmin = min(x(:)); xmax = max(x(:));  
ymin = min(y(:)); ymax = max(y(:)); zmin = min(z(:));

% 定义绘制圆锥体的位置
xrange = linspace(xmin,xmax,8);  yrange = linspace(ymin,ymax,8);  zrange = 3:4:15;  
[cx,cy,cz] = meshgrid(xrange,yrange,zrange);   

% 绘制圆锥体，并将缩放因子设置为5，使圆锥体大于默认大小。
figure
hcone = coneplot(x,y,z,u,v,w,cx,cy,cz,5);   

% 设置圆锥图颜色。
hcone.FaceColor = 'red';    
hcone.EdgeColor = 'none';

% 计算向量场的模(代表风速)，以生成用于切片命令的标量数据。
hold on
wind_speed = sqrt(u.^2 + v.^2 + w.^2);

hsurfaces = slice(x,y,z,wind_speed,[xmin,xmax],ymax,zmin);
set(hsurfaces,'FaceColor','interp','EdgeColor','none')
hold off

view(30,40);     daspect([2,2,1])	% 更改坐标区视图并设置数据纵横比

camlight right;     lighting gouraud	% 在相机的右侧添加光源，并设置光照
set(hsurfaces,'AmbientStrength',0.6);    hcone.DiffuseStrength = 0.8;

流线 streamline

clc;clear all;close all
load wind

% 标出流线起点
[startx,starty,startz] = meshgrid(80,20:10:50,0:5:15);
plot3(startx(:),starty(:),startz(:),'*r');

% 绘制流线
streamline(x,y,z,u,v,w,startx,starty,startz)
axis tight
view(3);

流管 streamtube

clc;clear all;close all

load wind
[sx,sy,sz] = meshgrid(80,[20 30 40],[5 10]);
verts = stream3(x,y,z,u,v,w,sx,sy,sz);
div = divergence(x,y,z,u,v,w);
streamtube(verts,x,y,z,-div);
view(3);
axis tight
shading interp
camlight 
lighting gouraud

流带 streamribbon

clc;clear all;close all
load wind
[sx,sy,sz] = meshgrid(80,[20 30 40],[5 10]);
verts = stream3(x,y,z,u,v,w,sx,sy,sz);
cav = curl(x,y,z,u,v,w);
spd = sqrt(u.^2 + v.^2 + w.^2).*.1;
streamribbon(verts,x,y,z,cav,spd);
axis tight
shading interp
view(3);
camlight; 
lighting gouraud

三维矢量场的综合表现

虽然我们说三维矢量场表示比较困难，但在MATLAB提供了如此多函数用以表示三维矢量场，根据要表示的场选择合适的函数，总能将自己想要表示的效果展现出来。