計算數學入門系列 - 計算區域的離散化（三）

情況在高維空間會變得有點複雜。在一維空間內， 我們上面考慮的都只是一段線段，沒有其他的可能性。在高維空間情況就有機會有點特殊。 在二維空間最簡單的情況是計算區域只是一個正方形或長方形。假設形狀是\([0,a]\times[0,b]\)，我們可以很簡單的把計算區域離散化

M=…;

N=…;

dx=a/M;

dy=b/N;

x=[0:dx:a];

y=[0:dy:b];

[x,y]=meshgrid(x,y);

可是如果計算區域並不是一個長方形，而是一個圓形，我們就不可能用一些\([0,dx]\times[0,dy]\) 的小長方形把這個形狀覆蓋。那到底我們應該如何把這個區域離散化 就有點複雜了。如果計算區域的形狀太過複雜 （比如說很多在工程上的應用），一般都會使用小三角形把整個計算區域覆蓋，我們把這個過程叫做三角剖分（Triangulation）。這些在有限元法（Finite Element Method） 的計算都是必須的。 如何把這些小三角形安放在整個計算區域內有時都不是那麼容易。有一些有限元法 的準確度會根據三角形的形狀而有所影響，我們可能希望三角形盡量是等邊三角形，也可能只希望這些小三角形的鈍角不會太細小。編寫這些程式有時候也不是太簡單。

這個三角剖分的過程，在三維空間還是可以做出來。可是高維度的計算，就非常複雜了。當然，就算計算區域是一些高維度的長方形，一般的偏微分方程計算還是非常困難，最主要的原因是這些取樣點的總數將會非常之多。如果每一方向我們有 N個點，在d維 空間，需要計算點的總數將會是\(N^d\)。 以d=3為例， 一般電腦都不可以計算到N=200 的情況。如果在四維空間，用相同的儲存空間， N就最多只可以是大約53。 對於大部份我們有興趣的情況，用53點根本不可能做出一個正確和精密的計算。在學術上，這個叫做維數災難 (Curse of Dimensionality)，一般我們做偏微分方程計算時都會盡量避免高於三維的計算。 

另外也有一些計算問題，在高位空間的計算區域並不是像上面那麼的直接。以一些動態界面問題為例子，我們希望計算氣泡在水裏面的運動情況。這個動態介面本身就是我們的計算區域。 最簡單的情況，是如何將取樣點平均分布在球體上，本身也已經是一個研究問題，叫做湯姆森問題（Thomson Problem） 。 有一些情況還是沒有找出來。 如果我們不需要平均分布，一般的三角剖分在這起泡表面上還是可以做出來。 但如果起泡在運動的過程形狀有所變化，在介面的一面黏合到他另外一邊，我們做計算時就必須要把這些在表面的三角形重新整合。在編寫程式是這些就非常困難。 所以就產生了另外一種方法，我們並不會把計算區域直接做離散化，而是把計算區域隱含在另外一個區域內。其中一個叫做水平值（Level Set）的方法，我們會定義一個函數 使到這個函數等於零的地方是我們需要的面。以一個球體為例，我們並不會直接在球體表面上取樣， 而是在這個球體存在的三維空間裏面作一個簡單的離散化。

N=…;

dx=3/(N-1);

x=[-1.5:dx:1.5];

[x,y,z]=meshgrid(x,x,x);

phi=sqrt(x.^2+y.^2+z.^2)-1;

p=patch(isosurface(x,y,z,phi,0));

isonormals(x,y,z,phi,p)

p.FaceColor = 'red';

p.EdgeColor = 'none';

daspect([1 1 1])

view(3)

camlight; lighting phong

在上面這段程式裏面，我們就定義了phi 這個函數，他等於零的地方就滿足了一個球體方程式 \(x^2+y^2+z^2=1\)。其他的指令，只是要求MATLAB 把這個函數等於零的地方找出來，然後把他畫成一個漂亮的圖。