利用matlab繪制二維均勻流線和向量場（向量場彩色箭頭，顏色隨變量變化）

• 0前言
• 1 均勻流線的繪制
• 2 繪制彩色的短線圖
• 3 繪制彩色的均勻流線
• 4 運動的彩色箭頭流線圖

0前言

之前一篇文章matlab流場可視化后處理我簡單介紹了很多繪制流場可視化的方法，但是在流線方面一直不是很滿意，因為matlab自帶的繪制流線函數(shù)streamslice沒有開源，有些參數(shù)（比如每根流線的顏色）不能修改等，所以我嘗試著自己寫一種繪制二維均勻流線的程序，以便滿足更多的需求。

其中，matlab官方file exchange里，也有一些大神的代碼給出了很好的結(jié)果：
Evenly Spaced Streamlines

Streamcolor

StreakArrow

關(guān)于如何用matlab繪制均勻的流線，我采用的方法是下面這篇論文提到的算法：
Jobard B., Lefer W. (1997) Creating Evenly-Spaced Streamlines of Arbitrary Density. In: Lefer W., Grave M. (eds) Visualization in Scientific Computing ’97. Eurographics. Springer, Vienna

1 均勻流線的繪制

論文Creating Evenly-Spaced Streamlines of Arbitrary Density中提出的繪制流線的算法可以描述為：

先初始化一根流線，然后在這個流線周圍逐漸的擴展新的流線。文章中定義了兩個距離，一個是用于產(chǎn)生新流線種子點的距離dstart，一個是用于結(jié)束流線生成的距離dend。新的流線擴展方式為，在距離舊的流線為dstart的隨機位置上，生成一個種子點，對該種子點進行新流線的生成（同時向頭和尾段生成流線），當(dāng)新流線兩端的點距某個舊流線的距離小于dend時，停止生成流線，將該流線放入舊流線內(nèi)保存。當(dāng)任意位置都不能滿足新流線的生成時，則認為所有流線都已經(jīng)生成，結(jié)束計算。

用偽代碼的形式可以表述為

1 隨機繪制第一根流線，作為初始流線。把初始流線加入已生成的流線列表內(nèi)。把初始流線記為當(dāng)前流線。 2 令Finish=False，當(dāng)Finish==False時，進行循環(huán)3 選擇一個種子點，該點距 當(dāng)前流線 的距離為dstart，且距離其它 已生成的流線列表內(nèi) 的流線距離大于dstart4 如果能夠找到一個符合條件的種子點，則生成新的流線。5 新流線以種子點為開始，向前向后伸展。當(dāng)流線的前端新計算的流線點距離某條流線的距離小于dend，則停止前端流線的生成。后端同樣。6 將新流線保存到已生成的流線列表內(nèi)。7 如果不能找到符合條件的種子點，則替換 當(dāng)前流線 ，重復(fù)48 如果 已生成的流線列表內(nèi) 所有流線都經(jīng)過步驟7嘗試一遍，則說明計算結(jié)束。令Finish==True， 停止2的循環(huán)。

由于在計算中，用到了大量的距離計算，比如生成新種子點時，需要不停的和其它所有流線計算距離，保證大于dstart，又比如在計算新流線時，新流線每次生成一個點，該點都需要和所有以有的流線進行距離計算，保證距離大于dend。

這種頻繁的距離計算數(shù)量，在量級上為N^2，這導(dǎo)致當(dāng)流線上點比較多或者流線比較多時，計算會非常緩慢。我們可以利用一些其它的方法減少這種運算，比如利用劃分網(wǎng)格的方式，使得每次計算距離，只需要計算周圍網(wǎng)格（比如2維的周圍網(wǎng)格有8個）內(nèi)所有點的距離，而無需計算全局所有點的距離。或者利用樹形搜索算法，可以將計算量降到NlgN的量級上。

當(dāng)然，這里我利用一種更簡單的方法進行計算。利用劃分網(wǎng)格的思想，將新種子點的生成位置限定在網(wǎng)格上。而且，假設(shè)如果兩個點在同一個網(wǎng)格內(nèi)，說明兩個點距離小于d，如果不在同一網(wǎng)格內(nèi)，兩個點距離大于d。這種假設(shè)雖然對計算結(jié)果不精準，但是仍然是可接受的，而且不需要計算距離。

用偽代碼的形式可以表述為

1 生成距離dstart的網(wǎng)格。生成距離dend的網(wǎng)格。網(wǎng)格內(nèi)初始化為0 2 當(dāng)start網(wǎng)格內(nèi)所有區(qū)域都為1，則停止循環(huán)3 尋找start網(wǎng)格區(qū)域為0的區(qū)域，隨機選擇一個區(qū)域作為生成種子點。4 以種子點為起始點，向前向后繪制流線5 查詢新生成的流線點對應(yīng)的end網(wǎng)格，記為當(dāng)前網(wǎng)格。6 如果當(dāng)前網(wǎng)格和上一個點的當(dāng)前網(wǎng)格是同一個，則認為當(dāng)前網(wǎng)格沒有移動，繼續(xù)生成下一個流線點7 如果當(dāng)前網(wǎng)格移動，且當(dāng)前end網(wǎng)格為1，則停止生成流線。8 如果當(dāng)前網(wǎng)格移動，且當(dāng)前end網(wǎng)格為0，則將當(dāng)前網(wǎng)格標記為1，繼續(xù)生成下一個流線點9 將流線保存。流線經(jīng)過的所有start網(wǎng)格標記為1。流線經(jīng)過的所有end網(wǎng)格標記為1。之后重復(fù)3。

matlab代碼為：

clear clc close all%載入數(shù)據(jù) load wind N=5; xx=x(:,:,N); yy=y(:,:,N); uu=u(:,:,N); vv=v(:,:,N); %由于這里R2018b之前的版本會出現(xiàn)兼容性問題，所以更改 %xmin=min(x,[],'all');xmax=max(x,[],'all'); %ymin=min(y,[],'all');ymax=max(y,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));%流線求解 streamline_sum=my_streamline(xx,yy,uu,vv,0.05);%繪制流線 figure(1) hold on xlim([xmin,xmax]); ylim([ymin,ymax]); xy_ratio = get(gca, 'DataAspectRatio'); xy_ratio = xy_ratio(1:2); xy_lim = axis;for j=1:size(streamline_sum,2)if size(streamline_sum{j},1)<=1continue%忽略異常流線endplot(streamline_sum{j}(:,1),streamline_sum{j}(:,2),'k');%繪制箭頭arrow_direction=streamline_sum{j}(end,:)-streamline_sum{j}(end-1,:);plot_arrow(xy_lim,xy_ratio,streamline_sum{j}(end,:),[0,0,0],0.8,arrow_direction) end hold offfunction streamline_sum=my_streamline(x,y,u,v,dstart) %0處理前設(shè)置 %設(shè)置網(wǎng)格密度(01區(qū)間內(nèi)歸一化的長度） %dstart=0.05;默認0.05 dend=0.5*dstart;%xmin=min(x,[],'all');xmax=max(x,[],'all'); %ymin=min(y,[],'all');ymax=max(y,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));%歸一化，將流場縮放為0-1區(qū)間的矩形 xn=(x-xmin)/(xmax-xmin); yn=(y-ymin)/(ymax-ymin); un=u/(xmax-xmin); vn=v/(ymax-ymin);num_start=ceil((0.5-dstart/2)/dstart)*2+1; num_end=ceil((0.5-dend/2)/dend)*2+1;%初始化所有網(wǎng)格點，0代表可以放置新點，1代表已經(jīng)存在原有的點 xy_start=zeros(num_start,num_start); xy_end=zeros(num_end,num_end);%1當(dāng)xy_start內(nèi)還有可放置的新點的位置時，進行循環(huán) k=0;%循環(huán)次數(shù)，也是流線個數(shù) while ~all(xy_start,'all')k=k+1;%2隨機一個start內(nèi)網(wǎng)格點作為種子點[start_id_y,start_id_x]=find(xy_start==0);randnum=randi(size(start_id_y,1));x_pos_i=id2axis(dstart,start_id_x(randnum,1));y_pos_i=id2axis(dstart,start_id_y(randnum,1));%3繪制流線streamline_i_1 = stream2(xn,yn, un, vn,x_pos_i,y_pos_i,0.2);streamline_i_2 = stream2(xn,yn,-un,-vn,x_pos_i,y_pos_i,0.2);%4以xy_end為標準，刪除自相交或間隔太近的點。并順便標記xy_end[streamline_i_1,xy_end,xy_start]=delete_self(streamline_i_1{1},xy_end,dend,xy_start,dstart);[streamline_i_2,xy_end,xy_start]=delete_self(streamline_i_2{1},xy_end,dend,xy_start,dstart);%5保存streamline_k=[flipud(streamline_i_2);streamline_i_1(2:end,:)];%新的流線streamline_sum{k}=[xmin+streamline_k(:,1)*(xmax-xmin),ymin+streamline_k(:,2)*(ymax-ymin)];%從歸一化還原 end endfunction xpoint=id2axis(distance,id) %取網(wǎng)格的中點 N=ceil((0.5-distance/2)/distance)*2+1;%分割的數(shù)量 min_distance=(1-(N-2)*distance)/2;%兩端最小的距離 if id==1xpoint=min_distance/2; elseif id==Nxpoint=1-min_distance/2; elsexpoint=min_distance+(id-1.5)*distance; end endfunction [sl_i,xy_end,xy_start]=delete_self(sl_i,xy_end,dend,xy_start,dstart) %sl_i streamline流線，兩列N行形式 N=size(sl_i,1); pos_id_last=axis2id(sl_i(1,1),sl_i(1,2),dend); xy_end(pos_id_last)=1;%第一個點標記%順便標記xy_start pos_id_s=axis2id(sl_i(1,1),sl_i(1,2),dstart); xy_start(pos_id_s)=1; for j=2:Npos_id_now=axis2id(sl_i(j,1),sl_i(j,2),dend);if pos_id_now~=pos_id_last%如果現(xiàn)在的點和原有的點在同一區(qū)域，則不管它%如果不在同一區(qū)域，檢測新的點是否已經(jīng)被占用if xy_end(pos_id_now)==1%如果該點被占用，說明出現(xiàn)與其它流線太近的情況，則直接停止j=j-1;breakelse%如果沒被占用，則把新點添加上xy_end(pos_id_now)=1;pos_id_last=pos_id_now;endend%順便標記xy_startpos_id_s=axis2id(sl_i(j,1),sl_i(j,2),dstart);xy_start(pos_id_s)=1; end sl_i(j:end,:)=[]; endfunction pos_id=axis2id(x,y,distance) N=ceil((0.5-distance/2)/distance)*2+1;%分割的數(shù)量 min_distance=(1-(N-2)*distance)/2;%兩端最小的距離 %x的位置 if x<=min_distancepos_id_x=1; elseif x>=1-min_distancepos_id_x=N; elsepos_id_x=ceil((x-min_distance)/distance)+1; end %y的位置 if y<=min_distancepos_id_y=1; elseif y>=1-min_distancepos_id_y=N; elsepos_id_y=ceil((y-min_distance)/distance)+1; end %xy轉(zhuǎn)ind pos_id=sub2ind([N,N],pos_id_y,pos_id_x); endfunction plot_arrow(xy_lim,xy_ratio,xy_arrow,arrow_color,arrow_width,arrow_direction) %初始化箭頭形狀（歸一化的形狀） arrow_0=[0,0;-1,0.5;-1,-0.5]; %對方向進行歸一化 a_dn=arrow_direction(:)./xy_ratio(:); a_dn=a_dn/sqrt(sum(a_dn.^2)); d=(xy_lim(4)-xy_lim(3)+xy_lim(2)-xy_lim(1))/2; %箭頭對窗口縮放 arrow_1=arrow_0*arrow_width*0.03*d; %箭頭旋轉(zhuǎn) arrow_2=arrow_1*[a_dn(1),a_dn(2);-a_dn(2),a_dn(1)]; %箭頭變形 xy_ratio_n=xy_ratio/sqrt(sum(xy_ratio.^2));%對比例尺歸一化 arrow_3=arrow_2.*xy_ratio_n+xy_arrow; fill(arrow_3(:,1),arrow_3(:,2),arrow_color,'EdgeColor','none') end

生成的流線圖為：

2 繪制彩色的短線圖

這個繪制思路很簡單，就是取每一個點，然后繪制短流線即可。
代碼如下：

clear clc close all %繪制短線彩色圖（長短相同）%載入數(shù)據(jù) load wind N=5; xx=x(:,:,N); yy=y(:,:,N); uu=u(:,:,N); vv=v(:,:,N); %xmin=min(xx,[],'all');xmax=max(xx,[],'all'); %ymin=min(yy,[],'all');ymax=max(yy,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));V2=sqrt(uu.^2+vv.^2);%計算速度%繪制變量顏色條 N_color=32; P=V2;%把速度作為變量 %P_max=max(P,[],'all');P_min=min(P,[],'all'); P_max=max(max(max(P)));P_min=min(min(min(P)));Num=numel(xx);%流線數(shù)量for k=1:Numstartx=xx(k);starty=yy(k);sl_i=stream2(xx,yy,uu,vv,startx,starty,[0.1, 20]); sl_sum{k}=sl_i{1};%保存流線 endmcp=colormap(parula(N_color));%生成顏色 [~,~,P_color] = histcounts(P(:),linspace(P_min,P_max,N_color));%顏色索引figure(1) hold on xlim([xmin,xmax]) ylim([ymin,ymax]) xy_ratio = get(gca, 'DataAspectRatio'); xy_ratio = xy_ratio(1:2); xy_lim = axis;for k=1:Numplot(sl_sum{k}(:,1),sl_sum{k}(:,2),'color',mcp(P_color(k),:))if size(sl_sum{k},1)<=1continueend%繪制箭頭arrow_direction=sl_sum{k}(end,:)-sl_sum{k}(end-1,:);plot_arrow(xy_lim,xy_ratio,sl_sum{k}(end,:),mcp(P_color(k),:),0.5,arrow_direction) end hold off caxis([P_min,P_max])%顏色條范圍 colorbarfunction plot_arrow(xy_lim,xy_ratio,xy_arrow,arrow_color,arrow_width,arrow_direction) %初始化箭頭形狀（歸一化的形狀） arrow_0=[0,0;-1,0.5;-1,-0.5]; %對方向進行歸一化 a_dn=arrow_direction(:)./xy_ratio(:); a_dn=a_dn/sqrt(sum(a_dn.^2)); d=(xy_lim(4)-xy_lim(3)+xy_lim(2)-xy_lim(1))/2; %箭頭對窗口縮放 arrow_1=arrow_0*arrow_width*0.03*d; %箭頭旋轉(zhuǎn) arrow_2=arrow_1*[a_dn(1),a_dn(2);-a_dn(2),a_dn(1)]; %箭頭變形 xy_ratio_n=xy_ratio/sqrt(sum(xy_ratio.^2));%對比例尺歸一化 arrow_3=arrow_2.*xy_ratio_n+xy_arrow; fill(arrow_3(:,1),arrow_3(:,2),arrow_color,'EdgeColor','none') end

繪制出的效果如下：

這些流線基本都是長短相同的，因為stream2函數(shù)生成的流線，流線上的點基本是等間距的。如果想做出速度V大的流線箭頭更長，速度小的流線更短的效果，這里用自己編寫的二階RK流線生成算法，進行流線的繪制。由于interp2速度太慢，這里用griddedInterpolant函數(shù)進行替換，優(yōu)化速度將近10倍（我也沒搞清楚為什么interp2速度這么慢）。

matlab代碼如下：

clear clc close all%繪制彩色短線圖（長短不同） %載入數(shù)據(jù) load wind N=5; xx=x(:,:,N); yy=y(:,:,N); uu=u(:,:,N); vv=v(:,:,N); %xmin=min(xx,[],'all');xmax=max(xx,[],'all'); %ymin=min(yy,[],'all');ymax=max(yy,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));V2=sqrt(uu.^2+vv.^2);%計算速度 %繪制變量顏色條 N_color=32; P=V2;%指定變量為V2 %P_max=max(P,[],'all');P_min=min(P,[],'all'); P_max=max(max(max(P)));P_min=min(min(min(P)));Num = numel(xx);num_streamline=20;%每條流線上的點數(shù)量 V2_max=max(max(max(V2))); dt=3.0*min([xx(1,2)-xx(1,1),yy(2,1)-yy(1,1)])/V2_max/num_streamline;for k=1:Numstartx=xx(k);starty=yy(k);sl_i=stream2_RK2(xx,yy,uu,vv,startx,starty,dt,num_streamline); sl_sum{k}=sl_i; endmcp=colormap(parula(N_color)); [~,~,P_color] = histcounts(P(:),linspace(P_min,P_max,N_color)); mlw=linspace(0.5,2,N_color);figure(1) hold on xlim([xmin,xmax]) ylim([ymin,ymax]) xy_ratio = get(gca, 'DataAspectRatio'); xy_ratio = xy_ratio(1:2); xy_lim = axis;for k=1:Numplot(sl_sum{k}(:,1),sl_sum{k}(:,2),'color',mcp(P_color(k),:),'linewidth',mlw(P_color(k)))if size(sl_sum{k},1)<=1continueend%繪制箭頭arrow_direction=sl_sum{k}(end,:)-sl_sum{k}(end-1,:);plot_arrow(xy_lim,xy_ratio,sl_sum{k}(end,:),mcp(P_color(k),:),0.5*mlw(P_color(k)),arrow_direction) end hold off caxis([P_min,P_max]) colorbarfunction streamline_i=stream2_RK2(x,y,u,v,startx,starty,dt,N) streamline_i=zeros(N,2); streamline_i(1,:)=[startx,starty]; x_old=startx;y_old=starty; F_u = griddedInterpolant(x',y',u','linear'); F_v = griddedInterpolant(x',y',v','linear'); for k=2:N%利用改進歐拉法（或者叫2階Runge-Kutta,預(yù)估校正）%interp2太慢，放棄 % u_K1=interp2(x,y,u,x_old,y_old,'linear')*dt; % v_K1=interp2(x,y,v,x_old,y_old,'linear')*dt;u_K1 = F_u(x_old,y_old)*dt;v_K1 = F_v(x_old,y_old)*dt; % u_K2=interp2(x,y,u,x_old+0.5*u_K1,y_old+0.5*v_K1,'linear')*dt; % v_K2=interp2(x,y,v,x_old+0.5*u_K1,y_old+0.5*v_K1,'linear')*dt;u_K2 = F_u(x_old+0.5*u_K1,y_old+0.5*v_K1)*dt;v_K2 = F_v(x_old+0.5*u_K1,y_old+0.5*v_K1)*dt;x_new=x_old+0.5*(u_K1+u_K2);y_new=y_old+0.5*(v_K1+v_K2);%保存streamline_i(k,:)=[x_new,y_new];x_old=x_new;y_old=y_new;if isnan(x_new) || isnan(y_new)streamline_i(k+1:end,:)=[];breakend end endfunction plot_arrow(xy_lim,xy_ratio,xy_arrow,arrow_color,arrow_width,arrow_direction) %初始化箭頭形狀（歸一化的形狀） arrow_0=[0,0;-1,0.5;-1,-0.5]; %對方向進行歸一化 a_dn=arrow_direction(:)./xy_ratio(:); a_dn=a_dn/sqrt(sum(a_dn.^2)); d=(xy_lim(4)-xy_lim(3)+xy_lim(2)-xy_lim(1))/2; %箭頭對窗口縮放 arrow_1=arrow_0*arrow_width*0.03*d; %箭頭旋轉(zhuǎn) arrow_2=arrow_1*[a_dn(1),a_dn(2);-a_dn(2),a_dn(1)]; %箭頭變形 xy_ratio_n=xy_ratio/sqrt(sum(xy_ratio.^2));%對比例尺歸一化 arrow_3=arrow_2.*xy_ratio_n+xy_arrow; fill(arrow_3(:,1),arrow_3(:,2),arrow_color,'EdgeColor','none') end

生成的流線效果如下：

3 繪制彩色的均勻流線

流線的繪制采用第1小節(jié)方法，但是采用控制stream2函數(shù)的點數(shù)來控制流線的長度（沒有采用第2小節(jié)的RK方法）。顏色的繪制采用第2小節(jié)的方法。

代碼如下：

clear clc close all%繪制彩色長線圖（長短不同） %載入數(shù)據(jù) load wind N=5; xx=x(:,:,N); yy=y(:,:,N); uu=u(:,:,N); vv=v(:,:,N); %xmin=min(xx,[],'all');xmax=max(xx,[],'all'); %ymin=min(yy,[],'all');ymax=max(yy,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));%流線求解 [streamline_sum,streamline_seed]=my_streamline_mutli(xx,yy,uu,vv,0.03,32);%繪制變量顏色條 V2=sqrt(uu.^2+vv.^2);N_color=16;%流線的顏色種類 P=V2;%把速度作為顏色的變量 %P_max=max(P,[],'all');P_min=min(P,[],'all'); P_max=max(max(max(P)));P_min=min(min(min(P)));mcp=colormap(jet(N_color)); P_seed=interp2(xx,yy,P,streamline_seed(:,1),streamline_seed(:,2));[~,~,P_color] = histcounts(P_seed(:),linspace(P_min,P_max,N_color)); mlw=linspace(0.5,2,N_color);%繪制流線 figure(1) hold on for j=1:size(streamline_sum,2)plot(streamline_sum{j}(:,1),streamline_sum{j}(:,2),'color',mcp(P_color(j),:)...,'linewidth',mlw(P_color(j)));%繪制流線 end xlim([xmin,xmax]) ylim([ymin,ymax]) caxis([P_min,P_max]) colorbar %繪制箭頭 xy_ratio = get(gca, 'DataAspectRatio'); xy_ratio = xy_ratio(1:2); xy_lim = axis; for k=1:size(streamline_sum,2)%繪制箭頭arrow_direction=streamline_sum{k}(end,:)-streamline_sum{k}(end-1,:);plot_arrow(xy_lim,xy_ratio,streamline_sum{k}(end,:),mcp(P_color(k),:),mlw(P_color(k)),arrow_direction) end hold offfunction [streamline_sum,streamline_seed]=my_streamline_mutli(x,y,u,v,dstart,num) %0處理前設(shè)置 %設(shè)置網(wǎng)格密度(01區(qū)間內(nèi)歸一化的長度） %dstart=0.05;默認0.05 dend=0.5*dstart;%xmin=min(x,[],'all');xmax=max(x,[],'all'); %ymin=min(y,[],'all');ymax=max(y,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));%歸一化，將流場縮放為0-1區(qū)間的矩形 xn=(x-xmin)/(xmax-xmin); yn=(y-ymin)/(ymax-ymin); un=u/(xmax-xmin); vn=v/(ymax-ymin);num_start=ceil((0.5-dstart/2)/dstart)*2+1; num_end=ceil((0.5-dend/2)/dend)*2+1;%初始化所有網(wǎng)格點，0代表可以放置新點，1代表已經(jīng)存在原有的點 xy_start=zeros(num_start,num_start); xy_end=zeros(num_end,num_end);%將流線劃分為num種，速度越大的流線越長 length_sl=linspace(5,40,num); V2=sqrt(un.^2+vn.^2); V2_max=max(max(max(V2))); V2_min=min(min(min(V2))); %V2_min=min(V2,[],'all');V2_max=max(V2,[],'all');V2_space=linspace(V2_min,V2_max,num+1);%1當(dāng)xy_start內(nèi)還有可放置的新點的位置時，進行循環(huán) k=0;%循環(huán)次數(shù)，也是流線個數(shù) while ~all(xy_start,'all')k=k+1;%2隨機一個start內(nèi)網(wǎng)格點作為種子點[start_id_y,start_id_x]=find(xy_start==0);randnum=randi(size(start_id_y,1));x_pos_i=id2axis(dstart,start_id_x(randnum,1));y_pos_i=id2axis(dstart,start_id_y(randnum,1));streamline_seed(k,:)=[x_pos_i,y_pos_i];%保存種子點V2_seed=interp2(xn,yn,V2,x_pos_i,y_pos_i);%計算種子點處的速度[~,~,sl_N] = histcounts(V2_seed,V2_space);num_streamline=round(length_sl(sl_N));%3繪制流線streamline_i_1 = stream2(xn,yn, un, vn,x_pos_i,y_pos_i,[0.1,num_streamline]);streamline_i_2 = stream2(xn,yn,-un,-vn,x_pos_i,y_pos_i,[0.1,num_streamline]);%4以xy_end為標準，刪除自相交或間隔太近的點。并順便標記xy_end[streamline_i_1,xy_end,xy_start]=delete_self(streamline_i_1{1},xy_end,dend,xy_start,dstart);[streamline_i_2,xy_end,xy_start]=delete_self(streamline_i_2{1},xy_end,dend,xy_start,dstart);%5保存streamline_k=[flipud(streamline_i_2);streamline_i_1(2:end,:)];%新的流線streamline_sum{k}=[xmin+streamline_k(:,1)*(xmax-xmin),ymin+streamline_k(:,2)*(ymax-ymin)];%從歸一化還原 end streamline_seed=[streamline_seed(:,1)*(xmax-xmin)+xmin,streamline_seed(:,2)*(ymax-ymin)+ymin]; endfunction xpoint=id2axis(distance,id) %取網(wǎng)格的中點N=ceil((0.5-distance/2)/distance)*2+1;%分割的數(shù)量 min_distance=(1-(N-2)*distance)/2;%兩端最小的距離 if id==1xpoint=min_distance/2; elseif id==Nxpoint=1-min_distance/2; elsexpoint=min_distance+(id-1.5)*distance; end endfunction [sl_i,xy_end,xy_start]=delete_self(sl_i,xy_end,dend,xy_start,dstart) %sl_i streamline流線，兩列N行形式 N=size(sl_i,1); pos_id_last=axis2id(sl_i(1,1),sl_i(1,2),dend); xy_end(pos_id_last)=1;%第一個點標記%順便標記xy_start pos_id_s=axis2id(sl_i(1,1),sl_i(1,2),dstart); xy_start(pos_id_s)=1; for j=2:Npos_id_now=axis2id(sl_i(j,1),sl_i(j,2),dend);if pos_id_now~=pos_id_last%如果現(xiàn)在的點和原有的點在同一區(qū)域，則不管它%如果不在同一區(qū)域，檢測新的點是否已經(jīng)被占用if xy_end(pos_id_now)==1%如果該點被占用，說明出現(xiàn)與其它流線太近的情況，則直接停止j=j-1;breakelse%如果沒被占用，則把新點添加上xy_end(pos_id_now)=1;pos_id_last=pos_id_now;endend%順便標記xy_startpos_id_s=axis2id(sl_i(j,1),sl_i(j,2),dstart);xy_start(pos_id_s)=1; end sl_i(j:end,:)=[]; endfunction pos_id=axis2id(x,y,distance) N=ceil((0.5-distance/2)/distance)*2+1;%分割的數(shù)量 min_distance=(1-(N-2)*distance)/2;%兩端最小的距離%x的位置 if x<=min_distancepos_id_x=1; elseif x>=1-min_distancepos_id_x=N; elsepos_id_x=ceil((x-min_distance)/distance)+1; end%y的位置 if y<=min_distancepos_id_y=1; elseif y>=1-min_distancepos_id_y=N; elsepos_id_y=ceil((y-min_distance)/distance)+1; end%xy轉(zhuǎn)ind pos_id=sub2ind([N,N],pos_id_y,pos_id_x); endfunction plot_arrow(xy_lim,xy_ratio,xy_arrow,arrow_color,arrow_width,arrow_direction) %初始化箭頭形狀（歸一化的形狀） arrow_0=[0,0;-1,0.5;-1,-0.5]; %對方向進行歸一化 a_dn=arrow_direction(:)./xy_ratio(:); a_dn=a_dn/sqrt(sum(a_dn.^2)); d=(xy_lim(4)-xy_lim(3)+xy_lim(2)-xy_lim(1))/2; %箭頭對窗口縮放 arrow_1=arrow_0*arrow_width*0.03*d; %箭頭旋轉(zhuǎn) arrow_2=arrow_1*[a_dn(1),a_dn(2);-a_dn(2),a_dn(1)]; %箭頭變形 xy_ratio_n=xy_ratio/sqrt(sum(xy_ratio.^2));%對比例尺歸一化 %arrow_3=arrow_2.*xy_ratio_n+xy_arrow;%由于兼容性問題，更改為下面語句 xy_ratio_n2=ones(3,1)*xy_ratio_n; arrow_3=arrow_2.*xy_ratio_n2+xy_arrow; fill(arrow_3(:,1),arrow_3(:,2),arrow_color,'EdgeColor','none') end

4 運動的彩色箭頭流線圖

一開始我本來是打算用第2小節(jié)的方式繪制一個動圖，但是發(fā)現(xiàn)原來的箭頭（粒子）離開之后，沒有新的箭頭進行補充。如果明確流場的入口和出口，可以采用出口離開多少粒子，就在入口添加多少粒子的方式，維持流場內(nèi)數(shù)量恒定。但是對于復(fù)雜流場則不容易做到這點。

因此，我采用了均勻流線的方式，每次運動完成后，刪除離得太近的流線，然后在空白處補充新的流線。效果如下：

這樣做的優(yōu)點在于每一幀都滿足均勻流線的要求，而且保證絕大多數(shù)流線不會憑空的產(chǎn)生和消失，即使是消失也是逐漸消失。

代碼如下：

clear clc close all %繪制彩色長線圖（長短不同，而且是會動的那種gif圖）%載入數(shù)據(jù) load wind N=5; xx=x(:,:,N); yy=y(:,:,N); uu=u(:,:,N); vv=v(:,:,N); %xmin=min(xx,[],'all');xmax=max(xx,[],'all'); %ymin=min(yy,[],'all');ymax=max(yy,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));dstart=0.03;%網(wǎng)格寬度%流線求解 [streamline_sum,streamline_seed]=my_streamline_mutli(xx,yy,uu,vv,dstart,32);%繪制變量顏色條 V2=sqrt(uu.^2+vv.^2);N_color=16;%流線的顏色種類 P=V2;%把速度作為顏色的變量 %P_max=max(P,[],'all');P_min=min(P,[],'all'); P_max=max(max(max(P)));P_min=min(min(min(P)));mcp=colormap(jet(N_color)); P_seed=interp2(xx,yy,P,streamline_seed(:,1),streamline_seed(:,2));[~,~,P_color] = histcounts(P_seed(:),linspace(P_min,P_max,N_color)); mlw=linspace(0.5,2,N_color);%繪制流線 figure(1) hold on for j=1:size(streamline_sum,2)plot(streamline_sum{j}(:,1),streamline_sum{j}(:,2),'color',mcp(P_color(j),:)...,'linewidth',mlw(P_color(j)));%繪制流線 end xlim([xmin,xmax]) ylim([ymin,ymax])%固定圖窗大小 caxis([P_min,P_max])%設(shè)定顏色條范圍 colorbar%顯示色條xy_ratio = get(gca, 'DataAspectRatio'); xy_ratio = xy_ratio(1:2); xy_lim = axis; for k=1:size(streamline_sum,2)%繪制箭頭if size(streamline_sum{k},1)<=2continueendarrow_direction=streamline_sum{k}(end,:)-streamline_sum{k}(end-1,:);plot_arrow(xy_lim,xy_ratio,streamline_sum{k}(end,:),mcp(P_color(k),:),mlw(P_color(k)),arrow_direction) end hold off%保存gif圖 frame=getframe(gca); im=frame2im(frame); [I,map]=rgb2ind(im,36); imwrite(I,map,'temp.gif','gif', 'Loopcount',inf,'DelayTime',0.08);%第一張%后續(xù)的變化 for k=1:60%延長dt秒之后的圖像streamline_sum=my_streamline_time(xx,yy,uu,vv,streamline_sum,dstart,0.05);figure(1)clfhold onfor j=1:size(streamline_sum,2)sl_i=streamline_sum{j};if isempty(sl_i)continueendL_sl_i=size(sl_i,1);P_sl=interp2(xx,yy,P,sl_i(round(L_sl_i/2),1),sl_i(round(L_sl_i/2),2));[~,~,P_color] = histcounts(P_sl,linspace(P_min,P_max,N_color));P_color(P_color==0)=1;plot(sl_i(:,1),sl_i(:,2),'color',mcp(P_color,:)...,'linewidth',mlw(P_color));%繪制流線%繪制箭頭if size(sl_i,1)<=2continueendarrow_direction=sl_i(end,:)-sl_i(end-1,:);plot_arrow(xy_lim,xy_ratio,sl_i(end,:),mcp(P_color,:),mlw(P_color),arrow_direction)endhold offxlim([xmin,xmax])ylim([ymin,ymax])%固定圖窗大小caxis([P_min,P_max])%設(shè)定顏色條范圍colorbar%顯示色條pause(0.5)%保存gif圖frame=getframe(gca);im=frame2im(frame);[I,map]=rgb2ind(im,36);imwrite(I,map,'temp.gif','gif','WriteMode','append','DelayTime',0.08); endfunction [streamline_sum,streamline_seed]=my_streamline_mutli(x,y,u,v,dstart,num) %繪制流線 %0處理前設(shè)置%設(shè)置網(wǎng)格密度(01區(qū)間內(nèi)歸一化的長度） %dstart=0.05;默認0.05 dend=0.5*dstart; %xmin=min(x,[],'all');xmax=max(x,[],'all'); %ymin=min(y,[],'all');ymax=max(y,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));%歸一化，將流場縮放為0-1區(qū)間的矩形 xn=(x-xmin)/(xmax-xmin); yn=(y-ymin)/(ymax-ymin); un=u/(xmax-xmin); vn=v/(ymax-ymin); F_u = griddedInterpolant(xn',yn',un','linear'); F_v = griddedInterpolant(xn',yn',vn','linear'); F_u_n = griddedInterpolant(xn',yn',-un','linear'); F_v_n = griddedInterpolant(xn',yn',-vn','linear'); num_start=ceil((0.5-dstart/2)/dstart)*2+1; num_end=ceil((0.5-dend/2)/dend)*2+1;%初始化所有網(wǎng)格點，0代表可以放置新點，1代表已經(jīng)存在原有的點 xy_start=zeros(num_start,num_start); xy_end=zeros(num_end,num_end); %將流線劃分為num種，速度越大的流線越長 length_sl=linspace(5,40,num); V2=sqrt(un.^2+vn.^2); %V2_min=min(V2,[],'all');V2_max=max(V2,[],'all'); V2_max=max(max(max(V2))); V2_min=min(min(min(V2)));V2_space=linspace(V2_min,V2_max,num+1);%1當(dāng)xy_start內(nèi)還有可放置的新點的位置時，進行循環(huán) k=0;%循環(huán)次數(shù)，也是流線個數(shù) while ~all(xy_start,'all')k=k+1;%2隨機一個start內(nèi)網(wǎng)格點作為種子點[start_id_y,start_id_x]=find(xy_start==0);randnum=randi(size(start_id_y,1));x_pos_i=id2axis(dstart,start_id_x(randnum,1));y_pos_i=id2axis(dstart,start_id_y(randnum,1));streamline_seed(k,:)=[x_pos_i,y_pos_i];%保存種子點%3繪制流線streamline_i_1=stream2_RK2(F_u ,F_v ,x_pos_i,y_pos_i,0.01,20);streamline_i_2=stream2_RK2(F_u_n,F_v_n,x_pos_i,y_pos_i,0.01,20);%4以xy_end為標準，刪除自相交或間隔太近的點。并順便標記xy_end[streamline_i_1,xy_end,xy_start]=delete_self(streamline_i_1,xy_end,dend,xy_start,dstart);[streamline_i_2,xy_end,xy_start]=delete_self(streamline_i_2,xy_end,dend,xy_start,dstart);%5保存streamline_k=[flipud(streamline_i_2);streamline_i_1(2:end,:)];%新的流線streamline_sum{k}=[xmin+streamline_k(:,1)*(xmax-xmin),ymin+streamline_k(:,2)*(ymax-ymin)];%從歸一化還原 end streamline_seed=[streamline_seed(:,1)*(xmax-xmin)+xmin,streamline_seed(:,2)*(ymax-ymin)+ymin]; endfunction xpoint=id2axis(distance,id) %取網(wǎng)格的中點 N=ceil((0.5-distance/2)/distance)*2+1;%分割的數(shù)量 min_distance=(1-(N-2)*distance)/2;%兩端最小的距離 if id==1xpoint=min_distance/2; elseif id==Nxpoint=1-min_distance/2; elsexpoint=min_distance+(id-1.5)*distance; end endfunction [sl_i,xy_end,xy_start]=delete_self(sl_i,xy_end,dend,xy_start,dstart) %sl_i streamline流線，兩列N行形式 N=size(sl_i,1); pos_id_last=axis2id(sl_i(1,1),sl_i(1,2),dend); xy_end(pos_id_last)=1;%第一個點標記%順便標記xy_start pos_id_s=axis2id(sl_i(1,1),sl_i(1,2),dstart); xy_start(pos_id_s)=1; for j=2:Npos_id_now=axis2id(sl_i(j,1),sl_i(j,2),dend);if pos_id_now~=pos_id_last%如果現(xiàn)在的點和原有的點在同一區(qū)域，則不管它%如果不在同一區(qū)域，檢測新的點是否已經(jīng)被占用if xy_end(pos_id_now)==1%如果該點被占用，說明出現(xiàn)與其它流線太近的情況，則直接停止j=j-1;breakelse%如果沒被占用，則把新點添加上xy_end(pos_id_now)=1;pos_id_last=pos_id_now;endend%順便標記xy_startpos_id_s=axis2id(sl_i(j,1),sl_i(j,2),dstart);xy_start(pos_id_s)=1; end sl_i(j:end,:)=[];%刪除停止之后剩余的流線 endfunction pos_id=axis2id(x,y,distance) N=ceil((0.5-distance/2)/distance)*2+1;%分割的數(shù)量 min_distance=(1-(N-2)*distance)/2;%兩端最小的距離 %x的位置 if x<=min_distancepos_id_x=1; elseif x>=1-min_distancepos_id_x=N; elsepos_id_x=ceil((x-min_distance)/distance)+1; end %y的位置 if y<=min_distancepos_id_y=1; elseif y>=1-min_distancepos_id_y=N; elsepos_id_y=ceil((y-min_distance)/distance)+1; end %xy轉(zhuǎn)ind pos_id=sub2ind([N,N],pos_id_y,pos_id_x); endfunction plot_arrow(xy_lim,xy_ratio,xy_arrow,arrow_color,arrow_width,arrow_direction) %初始化箭頭形狀（歸一化的形狀） arrow_0=[0,0;-1,0.5;-1,-0.5]; %對方向進行歸一化 a_dn=arrow_direction(:)./xy_ratio(:); a_dn=a_dn/sqrt(sum(a_dn.^2)); d=(xy_lim(4)-xy_lim(3)+xy_lim(2)-xy_lim(1))/2; %箭頭對窗口縮放 arrow_1=arrow_0*arrow_width*0.03*d; %箭頭旋轉(zhuǎn) arrow_2=arrow_1*[a_dn(1),a_dn(2);-a_dn(2),a_dn(1)]; %箭頭變形 xy_ratio_n=xy_ratio/sqrt(sum(xy_ratio.^2));%對比例尺歸一化 %arrow_3=arrow_2.*xy_ratio_n+xy_arrow;%由于兼容性問題，更改為下面語句 xy_ratio_n2=ones(3,1)*xy_ratio_n; arrow_3=arrow_2.*xy_ratio_n2+xy_arrow; fill(arrow_3(:,1),arrow_3(:,2),arrow_color,'EdgeColor','none') endfunction streamline_i=stream2_RK2(F_u,F_v,startx,starty,dt,N) %繪制流線，采用RK方法 %利用改進歐拉法（或者叫2階Runge-Kutta,預(yù)估校正） streamline_i=zeros(N,2); streamline_i(1,:)=[startx,starty]; x_old=startx;y_old=starty; for k=2:Nu_K1 = F_u(x_old,y_old)*dt;v_K1 = F_v(x_old,y_old)*dt;u_K2 = F_u(x_old+0.5*u_K1,y_old+0.5*v_K1)*dt;v_K2 = F_v(x_old+0.5*u_K1,y_old+0.5*v_K1)*dt;x_new=x_old+0.5*(u_K1+u_K2);y_new=y_old+0.5*(v_K1+v_K2);%保存streamline_i(k,:)=[x_new,y_new];x_old=x_new;y_old=y_new;if isnan(x_new) || isnan(y_new)streamline_i(k+1:end,:)=[];breakendif 0>x_new || 0>y_new || x_new>1 || y_new>1%由于插值結(jié)果都是在0-1區(qū)間內(nèi)的歸一化數(shù)據(jù)，所以超出邊界的刪除。streamline_i(k+1:end,:)=[];breakend end endfunction streamline_sum=my_streamline_time(x,y,u,v,streamline_sum,dstart,dt) %0處理前設(shè)置 %xmin=min(x,[],'all');xmax=max(x,[],'all'); %ymin=min(y,[],'all');ymax=max(y,[],'all'); xmin=min(min(min(x)));xmax=max(max(max(x))); ymin=min(min(min(y)));ymax=max(max(max(y)));%歸一化，將流場縮放為0-1區(qū)間的矩形 xn=(x-xmin)/(xmax-xmin); yn=(y-ymin)/(ymax-ymin); un=u/(xmax-xmin); vn=v/(ymax-ymin);streamline_sum(cellfun(@isempty,streamline_sum))=[];%刪除空數(shù)組 %排序，優(yōu)先保證長線段不被刪除 L_sl=zeros(length(streamline_sum),1); for k=1:length(streamline_sum)L_sl(k)=sum((streamline_sum{k}(end,:)-streamline_sum{k}(1,:)).^2); end [~,I]=sort(L_sl,'descend' ); streamline_sum(:,[1:length(streamline_sum)])=streamline_sum(:,[I]);%設(shè)置網(wǎng)格密度(01區(qū)間內(nèi)歸一化的長度） %dstart=0.05;默認0.05 dend=0.5*dstart;%每一條流線隨時間的變化 F_u = griddedInterpolant(xn',yn',un','linear'); F_v = griddedInterpolant(xn',yn',vn','linear'); F_u_n = griddedInterpolant(xn',yn',-un','linear'); F_v_n = griddedInterpolant(xn',yn',-vn','linear');num_start=ceil((0.5-dstart/2)/dstart)*2+1; num_end=ceil((0.5-dend/2)/dend)*2+1;for k=1:length(streamline_sum)sl_i=streamline_sum{k};if isempty(sl_i)continueendif any(isnan(sl_i(end,:)))sl_i(1,:)=[];elsesl_i(1,:)=[];if isempty(sl_i)continueendsl_i_n=(sl_i-[xmin,ymin])./[(xmax-xmin),(ymax-ymin)];%歸一化流線sl_i_now=stream2_RK2(F_u ,F_v ,sl_i_n(end,1),sl_i_n(end,2),dt,2);sl_i_n=[sl_i_n;sl_i_now(2,:)];sl_i=sl_i_n.*[(xmax-xmin),(ymax-ymin)]+[xmin,ymin];%反歸一化endstreamline_sum{k}=sl_i; end streamline_sum(cellfun(@isempty,streamline_sum))=[];%刪除空數(shù)組xy_start=zeros(num_start,num_start); xy_end=zeros(num_end,num_end); for k=1:length(streamline_sum)sl_i=streamline_sum{k};sl_i=flipud(sl_i);%刪除尾部，保留頭部sl_i_n=(sl_i-[xmin,ymin])./[(xmax-xmin),(ymax-ymin)];%歸一化流線%以xy_end為標準，刪除自相交或間隔太近的點。并順便標記xy_end[sl_i_n,xy_end,xy_start]=delete_self(sl_i_n,xy_end,dend,xy_start,dstart);sl_i=sl_i_n.*[(xmax-xmin),(ymax-ymin)]+[xmin,ymin];%反歸一化sl_i=flipud(sl_i);streamline_sum{k}=sl_i; end %之后還有一些區(qū)域太空，所以重新在那些區(qū)域生成流線 while ~all(xy_start,'all')k=k+1;%2隨機一個start內(nèi)網(wǎng)格點作為種子點[start_id_y,start_id_x]=find(xy_start==0);randnum=randi(size(start_id_y,1));x_pos_i=id2axis(dstart,start_id_x(randnum,1));y_pos_i=id2axis(dstart,start_id_y(randnum,1));%3繪制流線streamline_i_1=stream2_RK2(F_u ,F_v ,x_pos_i,y_pos_i,0.01,20);streamline_i_2=stream2_RK2(F_u_n,F_v_n,x_pos_i,y_pos_i,0.01,20);%4以xy_end為標準，刪除自相交或間隔太近的點。并順便標記xy_end[streamline_i_1,xy_end,xy_start]=delete_self(streamline_i_1,xy_end,dend,xy_start,dstart);[streamline_i_2,xy_end,xy_start]=delete_self(streamline_i_2,xy_end,dend,xy_start,dstart);%5保存sl_i_n=[flipud(streamline_i_2);streamline_i_1(2:end,:)];%新的流線sl_i=sl_i_n.*[(xmax-xmin),(ymax-ymin)]+[xmin,ymin];%反歸一化streamline_sum{k}=sl_i; endend

總結(jié)

以上是生活随笔為你收集整理的利用matlab绘制二维均匀流线和向量场的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。