非矩形區域上圖形繪制的MATLAB實現*

翟龍余

(江蘇財經職業技術學院基礎教學部，江蘇淮安223003)

MATLAB系統是目前世界上最流行的一種數值計算和圖形、圖像處理工具軟件，目前已廣泛應用于高等數學輔助教學中.比如，使用MATLAB軟件輔助教學，可以很容易繪制出復雜的立體圖形，把曲線、曲面的形成和變化過程準確地模擬出來，這對提高教學效率和培養學生的空間想象能力可起到事半功倍的效果.

如果要在空間直角坐標系下，繪制z=f(x，y)，其中(x，y) ∈［a，b］×［c，d］所代表的三維曲面圖形，可按如下步驟，用MATLAB進行繪圖［1］.

(1)生成平面網格坐標矩陣：先將［a，b］在x方向劃分出m個等分點，再將［c，d］在y方向劃分出n個等分點，最后利用meshgrid函數生成網格點的坐標，橫坐標放在x矩陣中，縱坐標放在y矩陣中，兩個矩陣都是n×m的矩陣；

(2)計算z坐標的矩陣：在上一步的基礎上利用MATLAB的數值計算功能計算出函數z=f(x，y)的函數值，生成z坐標的矩陣；

(3)繪制三維曲面：MATLAB提供了mesh函數和surf函數來繪制三維圖形，mesh函數用于繪制三維網格圖，surf函數用于繪制三維曲面圖.

高等數學中的函數z=f(x，y)的自然定義域不一定是矩形區域［a，b］×［c，d］，很多函數的自然定義域是一般形狀的平面區域，還有實際應用問題的研究往往不是在函數的自然定義域上研究，而是在一個與具體問題相關的非矩形區域上研究，這種非矩形區域上函數的三維曲面繪制在MATLAB的相關教材中都沒有介紹，本文就此結合數學中的典型實例研究了非矩形區域上空間三維曲面圖形的繪制.

1 非矩形區域上的MATLAB圖形繪制——自變量賦nan法

文獻［2］介紹了非矩形區域上MATLAB圖形繪制的一種方法：先把區域分割成若干個窄條矩形，利用MATLAB中繪制空間曲面的庫函數作出窄條矩形的圖形，最后再進行疊加，得到任意區域上的函數的圖形.但是，按照該法，如果區域的邊界線較為復雜時，所取點有可能落到區域外，如果該區域是函數的自然定義域時，還有可能繪制不出圖形.

在MATLAB中，有一個特殊的變量nan(Not a Number)，表示“非數”，MATLAB在處理數據時，會跳過數據中的nan，所有的MATLAB作圖函數都會忽略含有nan的數據點，實現一種所謂“鏤空”的效果［3］.受此啟發，可以通過如下方法作出邊界線已知的任意平面區域上的三維圖形：

(1)用上述矩形區域繪圖法步驟1的方法，生成一個矩形區域的網格坐標，要求該矩形區域包含所繪圖形的繪圖區域；

(2)剔除繪圖區域以外點：先計算出繪圖區域邊界線上函數值，然后通過MATLAB循環與條件判斷語句，一一判斷網格點是否在作圖區域內，如不在，把區域外的網格點的某個坐標值重新賦值為nan；

(3)按照上述矩形區域作圖法步驟2和3的方法進行繪圖，就可繪制出所需圖形.

上述繪圖方法關鍵之處是把繪圖區域外的點的一個坐標的值即函數的一個自變量的值重新賦值為nan，所以稱此繪圖方法為自變量賦nan法.

以x型區域為例，說明該繪圖法的編程過程：

由y=g(x)，y=h(x)，x=a，z=b圍成x型區域D(其中g(x)≥h(x))，現要求繪制二元函數z=f(x，y)，(x，y) ∈D的三維曲面圖，具體步驟如下：

(1)任取一矩形區域D0，使D0?D，如取D0=［a，b］×［gmax，hmin］，其中gmax是函數g(x) 的最大值，hmin是函數h(x)的最小值；

(2)針對D0區域，通過meshgrid命令生成網格點矩陣x和y；

(3)以網格點的橫坐標矩陣x，作為自變量，計算出相應的區域D的兩條邊界線上的函數值矩陣g和h；

(4)通過循環語句，比較矩陣y的元素與g和h相應的元素之間大小，若y(i，j) ＞g(i，j)，或y(i，j)＜h(i，j) 時，y(i，j) 重新賦值為 nan，即y(i，j)=nan.實際上此時縱坐標為y(i，j)所對應的點在區域D以外；

(5)計算函數z=f(x，y)的函數值；

(6)利用 mesh(x，y，z)(或 surf(x，y，z))命令繪制圖形，即可得到所需圖形.

實例1 繪制馬鞍面z=xy，(x，y)∈ {(x，y)||x|+|y|≤10}.

該圖在空間關于x=y面對稱，其繪圖區域在xοy面上既關于x軸對稱，又關于y軸對稱.下面的圖形繪制充分利用了其對稱性的特點.繪圖程序如下：

%選取矩形區域［0，10］×［0，10］生成網格點，該區域包含了繪圖區域的第一象限部分

x1=0：0.1：10；

y1=0：0.1：10；

n=length(x1)；

m=length(y1)；

［x1，y1］=meshgrid(x1，y1)；

%計算出第一象限邊界線上的函數值矩陣；

y=10-x1；

%以下循環將第一象限內繪圖區域外的點的縱坐標值重新賦值為nan.

for i=1：m；

for j=1：n；

if y1(i，j) ＞ y(i，j)；

y1(i，j)=nan；

end

end

end

%利用對稱性分別計算出第一、二、三、四象限區域上的函數值矩陣 z1、z2 、z3 、z4.

z1=x1.*y1；

z2=-x1.*y1；

z3=z1；

z4=z2；

%在同一圖形窗口繪制出一、二、三、四象限繪圖區域上的圖形，最終完成馬鞍面圖形的繪制.

mesh(x1，y1，z1)

hold on

mesh(-x1，y1，z2)

mesh(-x1，-y1，z3)

mesh(x1，-y1，z4)

得到如圖1所示馬鞍面，看到圖1就知道了它稱為馬鞍面的原因.

圖1 馬鞍面

圖2 半橢球面

如果函數z=f(x，y)的繪圖區域為y型區域，其繪圖方法與上類似，通過比較橫坐標的值，把繪圖區域外的點的橫坐標值重新賦值為nan即可.如果繪圖區域是更一般區域，可把其分割成若干個x型區域或y型區域，然后在同一圖形窗口分別用上述方法繪制出每個區域圖形即可.

2 函數z=f(x，y)的非矩形自然定義域上的三維曲面圖形繪制——函數值賦nan法

當函數z=f(x，y)的自然定義域為非矩形區域D時，顯然此時圖形繪制屬于非矩形區域上的圖形繪制，因而用上述自變量賦nan法即可繪制出其圖形.

但對于自然定義域為非矩形區域D的函數z=f(x，y)，多數情況下，如果給定一個矩形區域D0，D0包含函數z=f(x，y)自然定義域D，即D0?D.當(x，y)∈D0且(x，y)?D時，MATLAB軟件計算z=f(x，y)的函數值得到結果是虛數或inf，而MATLAB繪圖時，含inf坐標的數據點在圖中是不會顯示的(也不會出現錯誤信息)，如果此時我們用nan直接替換函數值矩陣中的虛數，再用mesh(x，y，z)(或surf(x，y，z))命令繪制圖形，可以更簡便地繪制出其自然定義域D上的圖形.該方法關鍵之處是把函數值矩陣中的虛數重新賦值為nan，所以稱此繪圖方法為函數值賦nan法，其具體步驟如下：

(1)任取一矩形區域D0，使D0包含自然定義域D；

(2)針對D0區域，通過meshgrid命令生成網格點矩陣x和y；

(3)直接計算函數z=f(x，y)的函數值；

(4)通過循環語句，一一判斷z(i，j)的虛部是否為，當其虛部不為時，z(i，j)重新賦值為nan，即z(i，j)=nan.實際上該步驟作用就是剔除自然定義域以外的點；

(5)利用 mesh(x，y，z)(或 surf(x，y，z))命令繪制圖形，即可得到所需圖形.

解 該函數的定義域為

y=-3：0.0025：3；

m=length(x)；% 計算x中元素的個數；

n=length(y)；

［x，y］=meshgrid(x，y)；% 生成矩形區域 的網格點坐標；

z=2.*(1-x.^2./4-y.^2./9).^(1/2)；% 計算函數值

for i=1：n；

for j=1：m；

if imag(z(i，j))～=0；% 查找z中的虛數；

z(i，j)=nan；% 將z中的虛數重新賦值為nan；

end

end

end

mesh(x，y，z)

其圖形如圖2所示為半橢球面.

注意點：以上兩種繪圖方法，由于受MATLAB計算精度以及計算機硬件條件的影響，在編程過程中x與y的步長選擇要適當，不要太小也不要太大.

［1］劉衛國.MATLAB程序設計與應用［M］.第二版.北京：高等教育出版社.2011.

［2］任明慧，張月蓮.MATLAB中矩形區域約束的突破［J］.湖南文理學院學報，2008，20(3)：9-11.

［3］鄭阿奇.MATLAB實用教程［M］.第二版.北京：電子工業出版社，2008.