Matlab編程避免使用循環語句的方法

李 冰，周德亮，邢澤寧

(遼寧師范大學 數學學院，遼寧 大連 116000)

Matlab編程避免使用循環語句的方法

李 冰，周德亮，邢澤寧

(遼寧師范大學 數學學院，遼寧 大連 116000)

充分利用Matlab豐富的矩陣運算功能，解決距離矩陣的生成、一次性生成N項和函數在多個節點的函數值、一次性繪制多個節點及多個幾何圖形的問題。同時進行不同編程方式在生成有限覆蓋子域圖形時所用時間的對比實驗，結果表明，這種矩陣化編程方法語言簡潔，且程序運行速度更快。

Matlab；矩陣化編程；矩陣生成；和函數計算；多圖形繪制

0 引言

一般編程語言在處理大量數據生成的數組或者矩陣問題時，多數采用循環語句實現[1]。Matlab語言不僅可以使用循環編程，同時可以利用其特殊的函數運算[2-3]，將矩陣作為變量的基本單位，充分利用矩陣化運算代替循環編程[4]。本文利用Matlab矩陣運算，結合研究中遇到的幾個問題，通過對距離矩陣的生成[5]、N項和函數的計算及有限覆蓋子域的繪制幾個實例，介紹了幾種解決這些實例的矩陣化編程方法。

1 距離矩陣生成

考慮如下形式矩陣的生成方法，在徑向基插值和配點法中經常用到距離矩陣，即：

傳統方法一般使用雙循環語句，而使用Matlab編寫程序可以避免循環。為了生成該矩陣，需要先構造如下形式的矩陣：

由于Ndgrid函數的功能是生成函數和插值所需的X和Y分量的網格數據，可利用Ndgrid函數生成矩陣Xc、Xr、Yc、Yr，如下：

x=0:0.1:1;y=0:0.1:1; [X,Y]=meshgrid(x,y);P=[X(:),Y(:)];

[Xc,Xr]=ndgrid(P(:,1));[Yc,Yr]=ndgrid(P(:,2));

A=sqrt((Xc-Xr).^2+(Yc-Yr).^2);

在區域[0,1]*[0,1]內布置節點x、y，其中X代表xi構成的向量，Y代表yi構成的向量，P代表將矩陣X和Y的元素按列排序。Ndgrid函數生成Xc、Xr、Yc、Yr矩陣后，利用Matlab對矩陣的“點”運算得到A，即Matlab將矩陣作為變量，對矩陣的每一個元素進行直接運算。

本文以生成2-范數距離矩陣為例。同理，當矩陣A中元素為任意范數形式時，可以求得任意N維空間的距離矩陣。

2 一類和函數計算

因此可編程如下，其中“%”部分為程序注釋。

N=50;rbf=inline('sqrt(1+(e*r).^2)','r','e') ;e= .006; %定義節點數及徑向基函數，e為參數。

x=linspace(0,5,N);%布置節點為列向量。

kxi=linspace(0,10,N);rand('state',4);c=rand(N,1); % 布置已知空間節點及系數c。

Pf=rbf(abs(repmat(x',1,N)-repmat(kxi,size(x'),1)),e)*c;

程序為：

N=50;rbf=inline('sqrt(1+(e*r).^2)','r','e');e=0.006;

kxi=linspace(0,10,N);rand('state',4);c=rand(N,1);

A=ones(size(x));E=eye(n); % 生成全1矩陣及單位矩陣。

m=10;n=20;x=fix(10*rand(m,n)); % 在[0,10]內生成m行n列的矩陣。

Pf=rbf(abs(repmat(x,1,N)-kron(kxi,A)),e)*kron(c,E);

其中，P(:,1)為Matlab編程語言，代表矩陣P的第一列，P(:,2)代表矩陣P的第二列。

其程序可編寫為：

rbf=inline('sqrt(1+(e*r).^2)','r','e') ;e= .006;

N=25;rand('state',4);P=fix(10*rand(N,2)); % 生成N行2列的矩陣P。

kxi=linspace(0,10,N);c=rand(N,1);

Pf=rbf(sqrt((kron(ones(1,N),P(:,1))-kron(kxi,ones(size(P(:,1))))).^2+...(kron(ones(1,N),P(:,2))-kron(kxi,ones(size(P(:,2))))).^2),ep)*c;

即Matlab矩陣化編寫程序如下：

N=50;rbf=inline('sqrt(1+(e*r).^2)','r','e') ;e= .006;

kxi=linspace(0,10,N);rand('state',4);c=rand(N,1);

x=0:2:10;y=5:5:20; [X,Y]=meshgrid(x,y) % 生成網格節點X,Y。

Pf=rbf(sqrt((repmat(X,1,N)-kron(kxi,ones(size(X)))).^2 ...+(repmat(Y,1,N)-kron(kxi,ones(size(Y)))).^2),ep) *kron(c,eye(length(x)))

Matlab函數在處理數據時采用矩陣運算，避免了對節點個數及空間維數的雙重循環編程。在徑向基函數配點法的實際應用中，這種網格節點編程方式經常用于散亂數據節點的布置。

3 網格、節點及子域邊界繪制

考慮如何不使用循環編程，一次性繪出如圖1和圖2所示的多個節點及多條直線圖形。

圖1 同時繪制多個節點

圖2 一次性繪制多條直線

利用Matlab編程可以實現一次性繪制多個幾何圖形。

x=0:0.1:1; y=0:0.1:1; [X,Y]=meshgrid(x,y);

figure(1); plot(X,Y,'bo','MarkerSize',5); figure(2); plot(X,Y,'k',Y,X,'k');

徑向基函數單位分解配點法的應用中，在布置多個節點后，需要選擇有限覆蓋子域。一般地，選擇矩形、圓、橢圓或球作為覆蓋圖形。

本文以橢圓覆蓋圖形為例，在[0,1]*[0,1]區域內一次性繪制多個以c為中心點的橢圓(見圖3)。

hold on

c=0.1:0.2:0.9; n=length(c); % n為橫坐標布置的中心點長度

[Xc,Yc]=meshgrid(c); XYc=[Xc(:),Yc(:)];

m=linspace(0,2*pi,n.^2)'; a=0.185; b=0.14; % m為極坐標下布置角度，a,b為參數

CX=kron(ones(1,n^2),Xc(:))'+kron(ones(1,n^2),a*cos(m)); %生成每個子域的邊界橫坐標

CY=kron(ones(1,n^2),Yc(:))'+kron(ones(1,n^2),b*sin(m)); %生成每個子域的邊界縱坐標

plot(CX,CY,'r',Xc,Yc,'k.');

圖3 橢圓形覆蓋域

圖4 時間對比

在繪制橢圓形覆蓋域時，使用Matlab中計時函數tic、toc對程序運行的時間進行檢測，同時利用如下循環程序進行對比實驗。時間對比如圖4所示。

for i=1:Nc

CX=XY(i,1)+a*cos(m); CY=XY(i,2)+b*sin(m);

plot(CX,CY,'r');

end

實驗對比表明，當中心點個數c增大，橢圓域個數增加時，Matlab將一組數據進行矩陣向量化編程，同時繪制多個幾何圖形，提高了程序的運行效率。

4 結語

本文用3個科研中遇到的應用實例問題，利用Matlab在矩陣數據處理上的優勢，實現了以矩陣化編程代替循環程序的過程，體現了Matlab矩陣化編程方法在提高程序運行效率上的優勢。

[1] 顧麗紅，李傳秀，吳少剛.矩陣乘法C語言程序設計案例探究[J].計算機教育，2016(1)：149-152.

[2] 楊亞輝.用Matlab深入學習和理解矩陣知識[J].現代電子技術，2007(6)：175-177.

[3] 鄢喜愛，楊金民，田華.Matlab在數據處理和繪圖中的應用[J].科學技術與工程，2006(6)：3631-3633.

[4] KANSA E J.Multiquadrics—a scattered data approximation scheme with applications to computational fluid dynamics II. Solutions to parabolic,hyperbolic and elliptic partial differential equations[J].Comput Math Appl,1990,19(8/9):147-161.

[5] 張志涌，楊祖櫻.MATLAB教程R2012a[M].北京：北京航空航天大學出版社，2010.

(責任編輯：孫 娟)

李冰(1990-)，女，遼寧昌圖人，遼寧師范大學數學學院碩士研究生，研究方向為科學計算可視化；周德亮(1960-)，男,遼寧沈陽人，博士，遼寧師范大學數學學院副教授、碩士生導師，研究方向為微分方程數值解法及應用；邢澤寧(1992-)，女，遼寧昌圖人，遼寧師范大學數學學院碩士研究生，研究方向為科學計算可視化。

10.11907/rjdk.162869

TP301

A

1672-7800(2017)003-0015-03