2維及3維空間上隱函數的等值面(線)圖
2維，3維的目的都一樣，就是做出隱函數表示的結構圖，將函數值為0的點視為表面並顯示出來，然后計算等值線(面)所圍之外的面積(體積)占整個空間的面積(體積)的百分比。
一.2維平面隱函數等值線
2維平面上的等值線圖使用contourf函數就可以實現。應該也有許多其他方法。
使用函數 \(z=\cos(x)*\cos(y)+0.5\)
1. 做出帶填充的等值(高)線圖
clear all;
clc;
% 給出定義域，生成網格。
x = 0:0.01:2*pi;
y = 0:0.01:2*pi;
[X, Y] = meshgrid(x, y);
% 給出隱函數表達式
Z = cos(X).*cos(Y)+0.5;
% 做等值線圖
ax = figure;
[M, C] = contourf(X, Y, Z);
axis off;
C.LineWidth = 1;
C.ShowText = 'on';
生成的圖像為：
根據二元函數求極值的方法，\(\frac{\partial f}{\partial x}\)，\(\frac{\partial f}{\partial y}\)都等於0，且\(\frac{\partial^2 f}{\partial x^2}\frac{\partial^2 f}{\partial y^2} - (\frac{\partial^2 f}{\partial x \partial y})^2 > 0\)，則x,y為函數極值點，若\(\frac{\partial^2 f}{\partial x^2}和\frac{\partial^2 f}{\partial y^2}\)都大於0，則為極小值，都小於0，則為極大值。
求出函數的極大值點分別為(0, 0)、\(\left(\pi,\pi\right)\)，函數值為1.5。極小值點分別為(0, \(\pi\))、(\(\pi\), 0)，函數值為-0.5。
從圖像上可以看出極大值與極小值點，且函數在x,y上的周期均為\(2\pi\)。
2. 做出函數值為0的帶填充等值(高)線
修改上面調用contourf的形式，順便把ShowText關掉。
...
[M, C] = contourf(X, Y, Z, [0, 0]);
...
C.ShowText = 'off';
得到的圖像為：
上面的圖像在邊界上並沒有做出等高線，因此，使用max函數為圖像增加邊界。
...
% 給出隱函數表達式
Z = max(max(cos(X).*cos(Y)+0.5, X.*(X-2*pi+0.01)), ...
Y.*(Y-2*pi+0.01));
...
之后做出的圖像為：
之所以使用\(X.*(X-2*pi+0.01)\)與\(Y.*(Y-2*pi+0.01)\)，是因為使用2pi時無法做出邊界，原因可能是matlab作圖時是按照類似的坐閉右開區間繪制的，因此把X,Y減小0.01。
還有一個問題是填充問題，我想實現的是等高線內部填充，而不是外部填充。
查了半天也沒找到辦法，通過colormap設置顏色梯度，里面永遠是白色，如下圖。
map = [0 0 0;
1 1 1];
colormap(ax, map);
得到的圖像為：
最后才發現，colormap填充的是等值線之間的圖像，因為我只畫了一條\(z=0\)的等值線，所以只能填充等值線以外的。
只要指定一個小於二元函數最小值的值以及所需的值，如下取[-0.5, 0]，使用colormap進行內部填充，-0.5為函數最小值。
修改代碼為：
...
[M, C] = contourf(X, Y, Z ,[-0.5, 0]); % 因為最小值為-0.5
axis off;
...
map = [0, 0.6, 1;
1, 1, 1];
colormap(ax, map);
得到圖像：
3. 得到等高線外部占全部面積的百分比
因為使用的離散數據，將大於0的所有節點求和，並除以所有節點即近似為空隙所占面積的百分比。
...
lenX = length(x);
lenY = length(y);
tot = lenX * lenY;
num_in_con = 0;
for i=1:1:lenX
for j=1:1:lenY
if (Z(i, j) > 0)
num_in_con = num_in_con + 1;
end
end
end
per = num_in_con / tot;
求得結果為per=0.8160;
二. 3維隱函數等值面
3維的隱函數等值面使用函數isosurface()、patch()以及isonormals()。
以三元函數\(f = (x^2+(\frac{9}{4}*y^2+\sqrt{z}^3-x^2*z^3-\frac{9}{80}*y^2*z^3\)為例，程序如下:
clear all;
clc
x = -10:0.05:10;
y = -10:0.05:10;
z = -10:0.05:10;
[X, Y, Z] = meshgrid(x, y, z);
f = (X.^2+(9/4)*Y.^2+Z.^2-1).^3-X.^2.*Z.^3-(9/80).*Y.^2.*Z.^3;
h = patch(isosurface(X, Y, Z, f, 0));
isonormals(X, Y, Z, f, h);
% 設置圖像屬性
h.FaceColor = 'red';
h.EdgeColor = 'none';
alpha(1);
view([1, 1, 1]);
axis equal;
axis off;
camlight('right');
lighting gouraud;
xlabel('x'); ylabel('y'); zlabel('z');
其中isosurface(X,Y,Z,f,0);獲取隱函數的等值面，返回包含面和頂點的結構體，這里的面指的是構成隱函數曲面的三角面片，而頂點即為三角形的三個頂點。
patch()由包含面和頂點的結構體繪制隱函數圖,返回h為Patch屬性，可以使用h改變圖像外觀。結果圖如下：
isonormals()基於函數值從新計算等值面的法向量，否則圖像不准確，不使用isonormals()函數的圖像如下：
