33402 從數學的角度看降臨堂

二、降臨堂屋頂的基本數學性質

$-\displaystyle\frac{4\sqrt{10}}{3}+2\ln\Big (\frac{\sqrt{82}-1}{\sqrt{82}+1}\Big ) +2\ln\Big (\frac{\sqrt{10}+3}{\sqrt{10}-3}\Big )$; 利用已知橫截面求體積法求得屋頂體積為 $12-\displaystyle\frac{4}{\sqrt{3}}$ (詳細計算過程請參考附註 1與附註 2)。

附註:

1. 由 (1) 式及微積分中常用的曲面表面積公式可知一面屋頂表面積 \begin{eqnarray*} &=& \iint\limits_R \sqrt{\Big (\frac{\partial z}{\partial x}\Big )^2 + \Big (\frac{\partial z}{\partial y}\Big )^2 + 1} ~dA \\ &=& \iint\limits_R \sqrt{\frac{1}{x^4} +1} ~dA \\ &=& \int_{\frac{1}{\sqrt{3}}}^3 \int_{-x}^x \frac{\sqrt{x^4+1}}{x^2} dy dx \\ &=& \int_{\frac{1}{\sqrt{3}}}^3 \bigg (\frac{\sqrt{x^4+1}}{x^2} \cdot 2x\bigg ) dx \\ &=& \int_{\frac{1}{3}}^9 \bigg (\frac{\sqrt{u^2+1}}{u}\bigg ) du ~\hbox{ (利用變數變換 }~ u = x^2) \\ &=& \int_{\tan^{-1}\frac{1}{3}}^{\tan^{-1}9} \frac{\sec^3\theta}{\tan\theta} d\theta ~\hbox{ (利用變數變換 }~ u = \tan\theta, ~0 \le \theta \le \frac{\pi}{2}) \qquad\qquad\qquad\quad \\ &=& \int_{\tan^{-1}\frac{1}{3}}^{\tan^{-1}9} \frac{1}{\cos^2\theta\cdot\sin\theta} d\theta \\ &=& \int_{\tan^{-1}\frac{1}{3}}^{\tan^{-1}9} \frac{\sin\theta}{\cos^2\theta\cdot\sin^2\theta} d\theta \\ &=& - \int_{\frac{3}{\sqrt{10}}}^{\frac{1}{\sqrt{82}}} \frac{1}{v^2\cdot(1-v^2)} dv ~\hbox{ (利用變數變換 }~ v = \cos\theta) \\ &=& - \int_{\frac{3}{\sqrt{10}}}^{\frac{1}{\sqrt{82}}} \bigg (\frac{1}{v^2} + \frac{1}{2(1+v)} + \frac{1}{2(1-v)}\bigg ) dv \\ &=& \bigg (\frac{1}{v} - \frac{1}{2}\ln|1+v| + \frac{1}{2}\ln|1-v|\bigg ) \bigg |_{\frac{3}{\sqrt{10}}}^{\frac{1}{\sqrt{82}}} \\ &=& \sqrt{82} - \frac{\sqrt{10}}{3} + \frac{1}{2} \ln \bigg (\frac{\sqrt{82}-1}{\sqrt{82}+1}\bigg ) + \frac{1}{2} \ln \bigg (\frac{\sqrt{10}+3}{\sqrt{10}-3}\bigg ). \end{eqnarray*} 據此, 四面屋頂表面積 $=4\sqrt{82}-\displaystyle\frac{4\sqrt{10}}{3}+2\ln\Big (\frac{\sqrt{82}-1}{\sqrt{82}+1}\Big ) +2\ln\big (\frac{\sqrt{10}+3}{\sqrt{10}-3}\Big )$。
2. 我們用已知橫截面求體積法求之。在高度為 $z$ 時, 橫截面為一正方形, 這個正方形邊長為 $\displaystyle\frac{2}{z}$, 因此截面面積為 $\displaystyle\frac{4}{z^2}$。 據此, 屋頂體積 $=\displaystyle\int_{\frac{1}{3}}^{\sqrt{3}} \displaystyle\frac{4}{z^2}dz =-\frac{4}{z}\Big |_{\frac{1}{3}}^{\sqrt{3}}=12-\frac{4}{\sqrt{3}}$。
3. 第四章第二節及 第六章第一節:
圓田術曰: 半周半徑相乘為積步。
徽注云: 按半周為從, 半徑為廣, 故廣從相乘為積步也。 假令圓徑二尺, 圓中容六觚之一面, 與圓徑之半, 其數均等。合徑率一而外周率三也。
又按: 為圖, 以六觚之一面乘一弧半徑, 三之, 得十二觚之冪。若又割之, 次以十二觚之一面乘一弧之半徑, 六之, 則得二十四觚之冪。 割之彌細, 所失彌少。割之又割, 以至於不可割, 則與圓周合體而無所失矣。觚面之外, 猶有餘徑, 以面乘餘徑, 則冪出弧表。 若夫觚之細者, 與圓合體, 則表無餘徑。表無餘徑, 則冪不外出矣。以一面乘半徑, 觚而裁之, 每輒自倍。 故以半周乘半徑而為圓冪。以此周、 徑, 謂至然之術, 非周三徑一之率也。$\cdots$

1. 方亭: 就是像圖七下方的方錐台。劉徽用了很漂亮的方法證明了方亭體積, 是讀者欣賞時不可錯過之處。 我們猜測當初這名稱應該是取亭子的屋頂上小下大的特色吧。
正方田: 正方形。
廣、 從: 即長方形的長、寬。
積、 積尺: 可能是面積或體積, 請參閱上下文。
答曰: 十丈: 這裡沒寫錯, 古代沒有平方單位、立方單位。
冪: 面積。
方堢壔: 正立方體。

1. 圖九是用 MATHEMATICA 畫的。指令如下:
Show[%,PlotRange $\rightarrow$ {1/3,1.732},Boxed $\rightarrow$ False,Axes $\rightarrow$ None]
2. 圖十是用 MATLAB 畫的。我們寫成指令檔 (commend file) 如下:
   clear
length=0.1;
[X,Y]=meshgrid(-3:length:3);
step=1+6/length;
for i=1:step
for j=1:step
if (abs(X(i,j))< 0.5) & (abs(Y(i,j))< 0.5)
Z(i,j)=0.0;
else
Z(i,j)=1.0/max(abs(X(i,j)),abs(Y(i,j)));
end
end
end
surf(X,Y,Z);
view([2.2 1.2 3]);
axis off

3. 圖十一是用 FORTRAN 求出斜邊的點, 再用 MATLAB 畫的:
FORTRAN 檔: (此檔可改用 MATLAB 檔, 寫起來會更加容易。)
    PROGRAM MAIN
IMPLICIT NONE
REAL,PARAMETER :: LENGTH=0.01
REAL,PARAMETER :: TOP=0.5
REAL,PARAMETER :: BOTTOM=3.0
REAL :: A,B,TOTLEN
INTEGER :: SIZE,I
TOTLEN=BOTTOM-TOP
SIZE=INT(TOTLEN/LENGTH)
OPEN(UNIT=100, FILE='l1.dat')
OPEN(UNIT=200, FILE='l2.dat')
OPEN(UNIT=300, FILE='l3.dat')
OPEN(UNIT=400, FILE='l4.dat')
DO I=0,SIZE
A=TOP+REAL(I)*LENGTH
B=1.0/A
WRITE(100,"(3(1X,F6.3))") -A,-A,B
WRITE(200,"(3(1X,F6.3))")  A,-A,B
WRITE(300,"(3(1X,F6.3))") -A, A,B
WRITE(400,"(3(1X,F6.3))")  A, A,B
END DO
STOP
END


MATLAB 檔:
      clear;
x1=[3 -3 -3 3 3];
y1=[3  3 -3 -3 3];
z1=[1/3 1/3 1/3 1/3 1/3];
x2=[0.5 -0.5 -0.5 0.5 0.5];
y2=[0.5  0.5 -0.5 -0.5 0.5];
z2=[1/0.5 1/0.5 1/0.5 1/0.5 1/0.5];
plot3(x1,y1,z1,'k-', ...
l1(:,1),l1(:,2),l1(:,3),'k-',l2(:,1),l2(:,2),l2(:,3),'k-', ...
l3(:,1),l3(:,2),l3(:,3),'k-',l4(:,1),l4(:,2),l4(:,3),'k-', ...
x2,y2,z2,'k-');
view([2.2 1.2 3]);
axis off

4. 圖十二也是用 FORTRAN 求出格子點, 再用 MATLAB 畫的:
FORTRAN 檔: (此檔可改用 MATLAB 檔, 寫起來會更加容易。)
    PROGRAM MAIN
IMPLICIT NONE
INTEGER,PARAMETER :: SIZE=600
REAL :: A,B,X,Y,LENGTH
REAL :: C(SIZE+1,SIZE+1)
INTEGER :: I,J,K
LENGTH=6.0/REAL(SIZE)
OPEN(UNIT=100, FILE='r.dat')
OPEN(UNIT=200, FILE='x.dat')
OPEN(UNIT=300, FILE='y.dat')
DO I=0,SIZE
DO J=0,SIZE
A=-3.0+REAL(I)*LENGTH
B=-3.0+REAL(J)*LENGTH
IF (ABS(A) > = 0.499 .OR. ABS(B) > = 0.499) THEN
C(I+1,J+1)= 1.0/MAX(ABS(A),ABS(B))
ELSE
C(I+1,J+1)= 0.0
END IF
WRITE(100,"(F6.3)") C(I+1,J+1)
END DO
END DO
DO K=0,SIZE
X=-3.0+REAL(K)*LENGTH
Y=-3.0+REAL(K)*LENGTH
WRITE(200,"(F6.3)") X
WRITE(300,"(F6.3)") Y
END DO
STOP
END


MATLAB 函數檔 (因為 FORTRAN 程式得到的結果必須用陣列存起來, 如果用矩陣形式存, 行太長會自動跳行。 因此, 我們要用函數檔將陣列改成方陣) :
    function A=arrtomat(x)
[k,l]=size(x);
m=sqrt(k);
for i=1:m
for j=1:m
A(i,j)=x(m*(i-1)+j,1);
end
end


MATLAB指令檔:
    clear