MATLAB繪制3D隱函數曲面的幾種方法

背景介紹
Matlab提供了一系列繪圖函數，常見的包括繪制2D曲線的plot函數、繪制2D隱函數曲線的ezplot函數、繪制3D曲面的mesh和surf函數、繪制3D顯函數曲面的ezmesh和ezsurf函數。值得注意的是，ez系列的繪圖函數里只有ezplot是繪制隱函數曲線的，ezmesh和ezsurf都是畫顯函數曲面的（不要被ez的名字誤解了）。遺憾的是，matlab里並沒有提供直接繪制3D隱函數曲面的函數。本帖的目的就是歸納總結幾種方便易用的繪制隱函數曲面的辦法。 

 

問題描述

如何繪制 3 元方程 f(x, y,z) = 0 確立的隱函數曲面 z = g(x,y) ？其中，方程 f(x, y,z) = 0 無法求解 z 關於 x 、 y 的表達式，即 g(x, y) 的顯式表達式無法獲取。 

 

准備工作——基礎函數介紹

為了解決上述問題，我們需要先 對幾個重要的圖形函數 isosurface 、 patch 、 isonormals 取得初步的了解，如果您已經對這三個函數很熟悉，可以直接跳過這一步。 

 

l.  isosurface 等值面函數

調用格式： fv = isosurface(X,Y,Z,V,isovalue)

作用：返回某個等值面（由 isovalue 指定）的表面（ faces ）和頂點（ vertices ）數據，存放在結構體 fv 中（ fv 由 vertices 、 faces 兩個域構成）。如果是畫隱函數 v = f(x,y,z) = 0 的三維圖形，那么等值面的數值為 isovalue = 0 。 

 

2.  patch函數

調用格式： patch(X,Y,C) 以平面坐標 (X, Y) 為頂點，構造平面多邊形， C 是 RGB 顏色向量

                    patch(X,Y,Z,C) 以空間 3-D 坐標 (X, Y,Z) 為頂點，構造空間 3D 曲面， C 是 RGB 顏色向量

                    patch(fv) 通過包含vertices、faces兩個域的結構體fv來構造3D曲面，fv可以直接由等值面函數isosurface得到 
例如：patch(isosurface(X,Y,Z,V,0))

 

3.  isonormals等值面法線函數

調用格式： isonormals(X,Y,Z,V,p)

實現功能：計算等值面 V 的頂點法線，將 patch 曲面 p 的法線設置為計算得到的法線（ p 是 patch 返回得到的句柄）。如果不設置法線的話，得到曲面在過渡地帶看起來可能不是很光滑 

 

有了上述三個函數后，我們已經具備間接繪制3D隱函數曲面的能力了。下面以方程
f(x,y, z) = x.*y.*z.*log(1+x.^2+y.^2+z.^2)-10 = 0為例，講解如何畫3D隱函數曲面。

 

解決辦法一：isosurface + patch+ isonormals

實現原理：先定義 3 元顯函數 v =f(x, y, z), 則 v = 0 定義的等值面就是 z = g(x,y) 的 3D 曲面。利用 isosurface 函數獲取 v= 0 的等值面，將得到的等值面直接輸入給 patch 函數，得出 patch 句柄 p ，並畫出 patch 曲面的平面視角圖形。對 p 用 isonormals 函數設置曲面頂點數據的法線，最后設置顏色、亮度、 3D 視角，得到 3D 曲面。

 

代碼如下：

 

f = @(x,y,z) x.*y.*z.*log(1+x.^2+y.^2+z.^2)-10; % 函數表達式 [x,y,z] = meshgrid(-10:.2:10,-10:.2:10,-10:.2:10); % 畫圖范圍 v = f(x,y,z); h = patch(isosurface(x,y,z,v,0)); isonormals(x,y,z,v,h) set(h,'FaceColor','r','EdgeColor','none'); xlabel('x');ylabel('y');zlabel('z'); alpha(1) grid on; view([1,1,1]); axis equal; camlight; lighting gouraud 

 

代碼說明：

• alpha函數用於設置patch曲面的透明度（可以是0~1任意數值），1 表示不透明，0 表示最大透明度。如果想設置透明度為0.7，可以修改alpha(1)為alpha(0.7)。

• 使用此代碼解決特定問題時，只需將第1行的函數表達式替換為特定問題的函數表達式，將第2行數據（x、y、z）范圍換成合適的范圍，后續代碼無需任何變動。

得到圖形： 

 

 

解決辦法二：Mupad

 

Mupad 符號引擎里提供了現成的三維隱函數畫圖函數： Implicit3d

 

在 matlab 里開啟 Mupad 的方法是：在 commandwindow 里輸入 mupad 來啟動一個 notebook 。在啟動的 notebook 里再輸入如下代碼：

 

 

plot(plot::Implicit3d(x*y*z*ln(1+x^2+y^2+z^2)-10, x = -10..10, y = -10..10, z = -10..10), Scaling = Constrained)

 

得到如下圖形：

 

解決辦法三：第三方工具包ezimplot3

 

在 matlab central 的 file exchange 上有一個非常優秀的繪制 3 維隱函數的繪圖函數，叫 ezimplot3 。感興趣的可以在如下鏈接下載：

http://www.mathworks.com/matlabcentral/fileexchange/23623-ezimplot3-implicit-3d-functions-plotter

 

 

ezimplot3 一共有三種參數調用方式： 

 

• ezimplot3(f) 畫函數f(X,Y,Z)= 0 在-2*pi< X < 2* pi, -2* pi < Y < 2* pi, -2* pi < Z < 2* pi上的圖形
• ezimplot3(f, [A,B])畫函數f(X,Y,Z)= 0 在A< X < B, A < Y < B, A < Z < B上的圖形
• ezimplot3(f, [XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX])畫函數f(X,Y,Z)= 0 在XMIN< X < XMAX, YMIN < Y < YMAX, ZMIN < Z < ZMAX上的圖形 
  
  

ezimplot3 使用方法：解壓 ezimplot3.zip ，將解壓得到的 ezimplot3.m 添加到 matlab 當前搜索路徑后就可以使用了。 代碼為：

f = @(x,y,z) x*y*z*log(1+x^2+y^2+z^2)-10; ezimplot3(f,[-10,10]); % [-10, 10] 表示圖形范圍x、y、z都在區間[-10, 10] 

 

若干說明： 

• ezimplot3和方法一本質上完全相同。即ezimplot3實際上也是基於isosurface+ patch + isonormals的實現
• ezimplot3與方法一的圖形視覺效果相同，唯一的區別是，ezimplot3的使用了0.7的透明度：alpha(0.7)
• ezimplot3在方法一基礎上增加了一些外包功能，如：允許函數句柄f是非向量化的函數（即函數定義無需.*  ./  .^），這在ezimplot3內部會自動調用vectorize實現函數向量化。另外，ezimplot3可以在調用的時候方便的設定坐標范圍。

 

常見問題和解決辦法：

 

• 常見問題：很多人在使用以上方法后，經常出現的問題是代碼沒有任何錯誤，程序可以運行，就是出來的圖形只有一個空坐標軸，看不到圖形。

• 問題分析：出現這種問題的原因是圖形的顯示區域沒設對。比如，我們上述三種方法都是在x為-10到10的范圍內，如果你設的范圍內本身就沒有圖形，那當然就看不到圖形了。
• 解決辦法：把圖形顯示范圍重新設置對即可，如果不知道圖形的大致范圍，就手工多改幾次，直到看到圖形為止

• 方法一，圖形范圍是在第2句的meshgrid函數決定的，meshgrid里給出的x、y、z范圍就是最終畫圖范圍，修改meshgrid語句即可。
• 方法二（Mupad），x =-10..10, y = -10..10, z = -10..10是表示顯示范圍，修改這里即可。
• 方法三，用ezimplot3(f,[A, B]) ezimplot3(f, [XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX])兩種方式控制圖形顯示范圍。

 

后記：slice切片函數

 

matlab 還提供一種畫切片圖形的函數 slice ， slice 做出的圖是在切片上用顏色表示 v 的值。有時，我們畫切片圖形也有助於我們理解一個 4 維圖形。 以   v= f(x,y,z) = x*y*z*exp(-(x^2+y^2+z^2))   為例，假設我們希望看 v =f(x,y,z) 在 x =0, y = 1, z = 1 這些平面切片的圖形，我們可以用以下代碼：

[x,y,z] = meshgrid(linspace(-2,2)); v = x.*y.*z.*exp(-(x.^2+y.^2+z.^2)); xslice = 0; yslice = 1; zslice = 1; slice(x,y,z,v,xslice,yslice,zslice) xlabel('x'); ylabel('y'); zlabel('z'); colormap hsv 

得到圖形為：