Matlab編程避免使用循環(huán)語句的方法

李 冰，周德亮，邢澤寧

(遼寧師范大學(xué) 數(shù)學(xué)學(xué)院，遼寧 大連 116000)

Matlab編程避免使用循環(huán)語句的方法

李 冰，周德亮，邢澤寧

(遼寧師范大學(xué) 數(shù)學(xué)學(xué)院，遼寧 大連 116000)

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

Matlab；矩陣化編程；矩陣生成；和函數(shù)計(jì)算；多圖形繪制

0 引言

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

1 距離矩陣生成

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

傳統(tǒng)方法一般使用雙循環(huán)語句，而使用Matlab編寫程序可以避免循環(huán)。為了生成該矩陣，需要先構(gòu)造如下形式的矩陣：

由于Ndgrid函數(shù)的功能是生成函數(shù)和插值所需的X和Y分量的網(wǎng)格數(shù)據(jù)，可利用Ndgrid函數(shù)生成矩陣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);

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

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

2 一類和函數(shù)計(jì)算

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

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

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

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

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

程序?yàn)椋?/p>

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]內(nèi)生成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) % 生成網(wǎng)格節(jié)點(diǎn)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函數(shù)在處理數(shù)據(jù)時(shí)采用矩陣運(yùn)算，避免了對節(jié)點(diǎn)個(gè)數(shù)及空間維數(shù)的雙重循環(huán)編程。在徑向基函數(shù)配點(diǎn)法的實(shí)際應(yīng)用中，這種網(wǎng)格節(jié)點(diǎn)編程方式經(jīng)常用于散亂數(shù)據(jù)節(jié)點(diǎn)的布置。

3 網(wǎng)格、節(jié)點(diǎn)及子域邊界繪制

考慮如何不使用循環(huán)編程，一次性繪出如圖1和圖2所示的多個(gè)節(jié)點(diǎn)及多條直線圖形。

圖1 同時(shí)繪制多個(gè)節(jié)點(diǎn)

圖2 一次性繪制多條直線

利用Matlab編程可以實(shí)現(xiàn)一次性繪制多個(gè)幾何圖形。

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');

徑向基函數(shù)單位分解配點(diǎn)法的應(yīng)用中，在布置多個(gè)節(jié)點(diǎn)后，需要選擇有限覆蓋子域。一般地，選擇矩形、圓、橢圓或球作為覆蓋圖形。

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

hold on

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

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

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

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

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

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

圖3 橢圓形覆蓋域

圖4 時(shí)間對比

在繪制橢圓形覆蓋域時(shí)，使用Matlab中計(jì)時(shí)函數(shù)tic、toc對程序運(yùn)行的時(shí)間進(jìn)行檢測，同時(shí)利用如下循環(huán)程序進(jìn)行對比實(shí)驗(yàn)。時(shí)間對比如圖4所示。

for i=1:Nc

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

plot(CX,CY,'r');

end

實(shí)驗(yàn)對比表明，當(dāng)中心點(diǎn)個(gè)數(shù)c增大，橢圓域個(gè)數(shù)增加時(shí)，Matlab將一組數(shù)據(jù)進(jìn)行矩陣向量化編程，同時(shí)繪制多個(gè)幾何圖形，提高了程序的運(yùn)行效率。

4 結(jié)語

本文用3個(gè)科研中遇到的應(yīng)用實(shí)例問題，利用Matlab在矩陣數(shù)據(jù)處理上的優(yōu)勢，實(shí)現(xiàn)了以矩陣化編程代替循環(huán)程序的過程，體現(xiàn)了Matlab矩陣化編程方法在提高程序運(yùn)行效率上的優(yōu)勢。

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

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

[3] 鄢喜愛，楊金民，田華.Matlab在數(shù)據(jù)處理和繪圖中的應(yīng)用[J].科學(xué)技術(shù)與工程，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].北京：北京航空航天大學(xué)出版社，2010.

(責(zé)任編輯：孫 娟)

李冰(1990-)，女，遼寧昌圖人，遼寧師范大學(xué)數(shù)學(xué)學(xué)院碩士研究生，研究方向?yàn)榭茖W(xué)計(jì)算可視化；周德亮(1960-)，男,遼寧沈陽人，博士，遼寧師范大學(xué)數(shù)學(xué)學(xué)院副教授、碩士生導(dǎo)師，研究方向?yàn)槲⒎址匠虜?shù)值解法及應(yīng)用；邢澤寧(1992-)，女，遼寧昌圖人，遼寧師范大學(xué)數(shù)學(xué)學(xué)院碩士研究生，研究方向?yàn)榭茖W(xué)計(jì)算可視化。

10.11907/rjdk.162869

TP301

A

1672-7800(2017)003-0015-03