next up previous
Next: 4 MATLAB 入門 Up: 2012年度情報処理II     第12回 MATLAB 体験 Previous: 2 レポート課題9 (4-b) について…ちょっと考えてみませんか

3 レポート課題9について

レポート課題9

(1)
変数 n に値が設定されているならば、 Cos[2Pi/n] とすると $ \cos\frac{2\pi}{n}$ が計算できるますが、 $ 20$ 個の値を計算するので、 Table[ ] を使うなどの工夫をすべきでしょう。
In[1] := Table[Cos[2Pi/n],{n,20}]
あるいは
In[2] := Table[{n,Cos[2Pi/n]},{n,20}]
5つずつ区切って表示すると

\begin{displaymath}
\begin{array}{ccccc}
1&-1&-\frac{1}{2}&0&\frac{1}{4} \left(...
...}{19}\right)&\sqrt{\frac{5}{8}+\frac{\sqrt{5}}{8}}.
\end{array}\end{displaymath}

知っている値が多いでしょう。簡単になっていませんが、 実は $ \cos\dfrac{\pi}{8}$ , $ \cos\dfrac{2\pi}{17}$ $ \sqrt{\mathstrut}$ を使って表せます。 $ \cos\dfrac{x}{2}>0$ であれば、

$\displaystyle \cos\frac{x}{2}=\sqrt{\frac{\cos x+1}{2}}
$

であるから、

$\displaystyle \cos\frac{\pi}{8}=\sqrt{\frac{\cos\frac{\pi}{4}+1}{2}}
=\frac{\sqrt{2+\sqrt{2}}}{2}.
$

これを確かめるには、引き算して FullSimplify[ ] すると良い:
In[3] := FullSimplify[Cos[Pi/8]-Sqrt[2+Sqrt[2]]/2]
注意すべきは FullSimplify[Cos[Pi/8]] としても、 $ \dfrac{\sqrt{2+\sqrt{2}}}{2}$ は得られないことである。 一方、 $ \cos\frac{2\pi}{17}$ ですが、実は

$\displaystyle \cos\frac{2\pi}{17}
=\frac{1}{16}
\left(
-1+\sqrt{17}+\sqrt{34...
...}}
+2\sqrt{17+3\sqrt{17}-\sqrt{34-2\sqrt{17}}-2\sqrt{34+2\sqrt{17}}}
\right)
$

です1
In[4] := a = (-1 + Sqrt[17] + Sqrt[34 - 2 Sqrt[17]] + 
     2 Sqrt[17 + 3 Sqrt[17] - Sqrt[34 - 2 Sqrt[17]] - 
        2 Sqrt[34 + 2 Sqrt[17]]])/16 - Cos[2 Pi/17]
In[5] := FullSimplify[a]
としてみて下さい。
(2)
計算すべきものはみな $ \dsp\sum_{k=1}^n\frac{1}{2^k}$ という形なので、 それを計算する関数を定義するのが第一歩。
In[6] := mysum[n_]:=Sum[1/2^k,{k,1,n}]
以下 mysum[3], mysum[5], mysum[10], mysum[50] と4回コマンドを入力しても良いですが、 リストを使って
In[7] := a=mysum[{3,5,10,50}]

$\displaystyle \frac{7}{8},
\frac{31}{32},
\frac{1023}{1024},
\frac{1125899906842623}{1125899906842624}
$

小数点以下 $ 60$ 桁表示すると
In[8] := N[a,60+1]

$\displaystyle \begin{array}{l}
0.8750000000000000000000000000000000000000000000...
...\
0.9999999999999991118215802998747676610946655273437500000000000.
\end{array}$

$ n$ 項までの和は、小数点以下 $ n$ 桁目まであれば十分で、 $ n+1$ 桁目より下は 0 のようですね。 ちょっと考えると理由が分かると思います。$ 10^n$ をかけてみると、
In[8]:= mysum10n[n_] := 10^n Sum[1/2^k, {k, 1, n}]
In[9]:= mysum10n[{3,5,10,50}]

$\displaystyle 875, 96875, 9990234375, 99999999999999911182158029987476766109466552734375.
$

つまり

$\displaystyle \frac{875}{10^3}, \frac{96875}{10^5},
\frac{9990234375}{10^{10}},
\frac{99999999999999911182158029987476766109466552734375}{10^{50}}.
$

(3)
In[] := a=3
In[] := x[1]=1
In[] := x[n_]:=x[n]=(x[n-1]+a/x[n-1])/2
In[] := tablex=Table[x[n],{n,1,10}]
In[] := N[tablex,50+1]
tablex はとんでもないことになっているので (分母と分子がどでかい有理数)、 最初から小数計算すべきかもしれません。 x[1]=1 のかわりに x[1] =N[1,50+1] として始めるという手があります。 いずれにしても

$\displaystyle \begin{array}{l}
1.0000000000000000000000000000000000000000000000...
...525381038, \\
1.73205080756887729352744634150587236694280525381038
\end{array}$

第8項目までで、小数点以下$ 50$ 桁まで求まっていそうです。 100桁まで計算して、誤差をカンニングで調べてみましょう。
In[] := N[Abs[N[tablex,100]-Sqrt[3]]]

$\displaystyle 0.732051, 0.267949, 0.0179492, 0.0000920496,
2.44585\times10^{-9}, 1.72691\times 10^{-18},
8.6089\times10^{-37}, 2.13946\times10^{-73}, 0, 0
$

となります。大雑把に言って、 合致する桁数 (1, 2, 4, 9, 18, 37, 73 と) が倍々ゲームで増えていくことが分かります。 誤差の常用対数をプロットして見ると
\includegraphics[width=10cm]{eps/errornewton.eps}
(4)
(a)
楕円面 $ \dfrac{x^2}{2}+\dfrac{y^2}{3}+\dfrac{z^2}{4}=1$ を歪めずに 描くには BoxRatios->Automatic とするのが簡単でしょう。

一度に描いてしまうことも出来ます。 みんな等値面 (レベルセット) として描いてしまう:
In[] := ContourPlot3D[{x^2/2 + y^2/3 + z^2/4 == 1,
                      x + y + z == 3, x + y + z == -3},
                      {x, -2, 2}, {y, -2, 2}, {z, -7, 7}, 
                      BoxRatios -> Automatic]
あるいは、みんな関数のグラフとして描いてしまう:
In[] := Plot3D[{Sqrt[4 (1 - x^2/2 - y^2/3)],
               -Sqrt[4 (1 - x^2/2 - y^2/3)], -x - y -  3, -x - y + 3},
               {x, -2, 2}, {y, -2, 2}, BoxRatios -> Automatic]
しかし、別々に描いて合成する方法をマスターすることを勧めたいです (考え方として自然だと思います)。 別々に描いてよいならば、曲面ごとに描く方法を選択できて、 自由度が大きくなるので。

     平面はグラフとして、楕円面は等値面として描く:
In[] := L=2;
In[] := g1=Plot3D[{-x-y-3,-x-y+3},{x,-L,L},{y,-L,L},BoxRatios->Automatic]
In[] := g2=ContourPlot3D[x^2/2+y^2/3+z^2/4==1,{x,-L,L},{y,-L,L},{z,-L,L}]
In[] := g=Show[g1,g2]
この場合、Show[ ] の引数の順番は大事です。 大きいものを最初に持ってきます。

     しかし、PlotRange->All というオプションに気がつけば、 細かいことを気にする必要はありません。 一つ一つ描いて、最後に順番を気にしないでまとめます:
In[] := g1=Plot3D[-x-y-3,{x,-L,L},{y,-L,L}]
In[] := g2=Plot3D[-x-y+3,{x,-L,L},{y,-L,L}]
In[] := g3=ContourPlot3D[x^2/2+y^2/3+z^2/4==1,{x,-L,L},{y,-L,L},{z,-L,L}]     
In[] := Show[g1,g2,g3,PlotRange->All, BoxRatios->Automatic]

(b)
軸が $ z$ 軸で、頂点が原点である直円錐を描くことにします。 $ z=a r$ ($ a$ は定数) から、グラフとしての表現

$\displaystyle z=a\sqrt{x^2+y^2},
$

パラメーター曲面としての表現、

$\displaystyle (x,y,z)=(r\cos\theta,r\sin\theta,a r),
$

方程式の解集合 (≒関数のレベルセット) としての表現

$\displaystyle x^2+y^2=\left(\frac{z}{a}\right)^2,$   あるいは$\displaystyle \quad
x^2+y^2-\left(\frac{z}{a}\right)^2=0.
$

前のページに描いたように、 $ x=f(z)$ $ z$ 軸のまわりに回転してできる曲面という考え方をすると、 方程式は $ x^2+y^2=f(z)^2$ となります。 $ f(z)=\dfrac{z}{a}$ とすると、 上の方程式が得られます。
g = Plot3D[3 Sqrt[x^2 + y^2] Boole[x^2 + y^2 <= 1],
          {x, -1, 1}, {y, -1, 1}, BoxRatios -> Automatic]
g = ParametricPlot3D[{r Cos[t], r Sin[t], 3 r}, {r, 0, 1}, {t, 0, 2 Pi}]
g = ContourPlot3D[x^2 + y^2 - (z/3)^2 == 0,
                 {x, -1, 1}, {y, -1, 1}, {z, -3, 3}, BoxRatios -> Automatic]


next up previous
Next: 4 MATLAB 入門 Up: 2012年度情報処理II     第12回 MATLAB 体験 Previous: 2 レポート課題9 (4-b) について…ちょっと考えてみませんか
桂田 祐史
2012-07-11