Mathematica程式——計算函數對應映像所圍成的面積

 【問題】試求如下曲線所圍成圖形的面積。

x^2 + y^2 ≤ 16; x^2/36 + y^2  ≤ 1; (x - 2)^2 + (y+1) ≤ 9;

第一種演算法：

Area[x1_, x2_, y1_, y2_, step_] :=Block[{nx, ny, i, j, xx, yy, s = 0},    nx = IntegerPart[(x2 - x1)/step];    ny = IntegerPart[(y2 - y1)/step];    For[i = 1, i ≤ nx, i++,      For[j = 1, j ≤ ny, j++, xx = x1 + i*step; yy = y1 + j*step;        If[xx^2 + yy^2 ≤ 16 && xx^2/36 + yy^2 ≤ 1 && (xx -     2)^2 + (yy + 1)^2 ≤ 9, s = s + step^2]]];    Return[s]]

【運行結果】

Area[-2, 6, -2, 2, 0.01]8.8312Area[-1, 4, -1, 1, 0.001]8.83922

第一種方法的思想很簡單，就是細分矩形地區，將落在目的地區域內的小矩形的面積進行累加

面積的近似程度與細分的精度step有關，演算法較為笨拙，計算時間長，演算法複雜度為O(1/step^2)。

第二種方法：基於機率思想的蒙特卡羅演算法的啟迪

AreaMC[x1_, x2_, y1_, y2_, n_] :=Block[{rx, ry, i, s, calin = 0},    For[i = 1, i ≤ n, i++,      rx = Random[Real, {x1, x2}];      ry = Random[Real, {y1, y2}];      If[rx^2 + ry^2 ≤ 16 && rx^2/36 +       ry^2 ≤ 1 && (rx - 2)^2 + (ry + 1)^2 ≤ 9, calin = calin + 1]];    s = (x2 - x1)*(y2 - y1)*N[calin/n];    Return[s]]

 

AreaMC[-2, 6, -2, 2, 1000000]8.85734

 

該演算法是近似演算法，演算法的複雜度依賴於隨機產生的點的個數n，複雜度為O(n)