Chapter 15
章末「研究」のヒントと略解

各章の章末にある「研究」へのヒント，略解，あるいはコメントです．

注意． 本章は「補講」，すなわちオマケです．13 章までの内容は既知とさせていただ
きます．説明も親切さに欠ける部分もあるかもしれませんがご了承ください．

2.6 研究(1) 略．
√√(2) a + b の形の二重根号は明らかに外せそうな場合でも a, b が整数のときは√√良いが，有理数になるとだめになったりする．たとえば3 + 2 2 は外せる√√2が 3+2はダメ，といった具合である．Mathematica の FullSimplify が10
どのように式の簡略化をアルゴリズム化しているのか，実際のところはよくわからない．

(3) ちなみに (1 − i)i は ei log(1−i) で定義され，さらに log(1 − i) は ez = 1 − i を満たす無限個の複素数のことである．したがって (1 − i)i も無限個の複素数となるが，N[(1-I)^I] としてもひとつの数値しか出てこない．
私の専門は複素解析なので，こういう複素数の指数・対数にともなう「多価性」の問題はいつも気になるところ．数値的に確かめたりドキュメントセンターを見たりして，間違った結果にならないよう注意を払っている．*1*1 ちなみに指数関数 ez = exp(z) は無限級数で定義される一価関数なので問題ない．ただし ez を「e = 2.718 · · · の複素数 z 乗」と解釈すると，無限個の複素数が出てきて注意が必要である．

215

章末「研究」のヒントと略解

3.5 研究(1) Table[i + j, {i, 1, m}, {j, i, n}}](2) Table[i + j + k, {i, 1, n}, {j, i, n}, {k, 1, n}}](3) Table[{n, Floor[1 + Log[10, Sum[4 8^(i-1), {i, 1, n}]]]},{n, 1, 100}]//TableFormと入力してみよう．n と Sn の桁数の組が n = 1 から 100 まで表示されるはずである．S99 は 90 桁，S100 は 91 桁であり，1 ≤ d ≤ 91 を満たす自然数についてちょうど d 桁となる Sn が存在することが確認できる．では練習問題：
一般に，任意の自然数 d について，ちょうど d 桁となる Sn が存在するか？

(4) 略．

4.7 研究
Reduce は代数的な式変形で論理的に方程式を解いてくれる．厳密だが数値解析には向かない．FindRoot にはニュートン法が組み込まれている．複素数解を求めることもできるが，いつでも解が求まるとは限らない．（14 章参照．）

5.6 研究(1) 1 次式なので f (x) = ax + b とおく．さらに 左辺 − 右辺の関数を g(x) とおけば，g(0) = g(1) = g(2) = 0 が成り立つはずである．よって次のようにすれば a, b, c が求まる：

f[x_] = a x + b;g[x_] = (Integrate[f[y], {y, 0, x}]+ Integrate[(x + y)^2*f[y], {y, 0, 1}]) - (x^2 + c);
Solve[{g[0] == 0, g[1] == 0, g[2] == 0}, {a, b, c}]
答えは a = 3/2，b = −1/2, c = 5/24．SolveAlways を用いてSolveAlways[ Integrate[f[y], {y, 0, x}]+ Integrate[(x + y)^2*f[y], {y, 0, 1}]== x^2 + c, x]としても同じ結果が得られる．

(2) f[x_] = a x^3 + b x^2 + c x + d;3

sol = Solve[{f[1] == 1, f[-1] == -1,Integrate[b x^2 + c x + d, {x, -1, 1}] == 1},{b, c, d}]とすれば b = −3/4, c = 1 − a, d = 3/4 を得る．この結果を f に代入した式（g とおく）を考えて，

g = f[x]/.sol[[1]]gxx = D[g, x, x]
h = Integrate[gxx^2, {x, -1, 1/2}]とすれば，積分 I （ここでは h とおいた）は 27/8 + (27a)/4 + (27a2 )/2 となることがわかる．これは a の 2 次関数（下に凸）であるから，a に関する微分が 0 になるところで最小である．よってSolve[D[h, a] == 0, a]を実行して，a = −1/4 を得る．
(3) (a) は (2) と同様に解く：f[x_] = a x^2 + b x + c;h[x_] =(Integrate[f[t], {t, 0, x}]+ x Integrate[f[t], {t, -1, 1}]- (1/3)*(f[1] - f[-1]) )- (4 x^3 + p x^2 - 10 x - 4);
sol = Solve[{h[0] == 0, h[1] == 0, h[2] == 0, h[3] == 0},{a, b, c, p}]
これより a = 12, b = 6, c = −6, p = 3 を得る．SolveAlways を用いて
sol = SolveAlways[Integrate[f[t], {t, 0, x}]+ x * Integrate[f[t], {t, -1, 1}]- 1/3(f[1] - f[-1])== 4 x^3 + p x^2 - 10 x - 4, x]
としてもよい．(b) は問題 5.1 と同様に解く．
g = (x + 1) f[x] /.sol[[1]] とおくと，g(x) = (1 + x)(−6 + 6x + 12x2 )となっていることがわかる．極値を求めるために微分して 0 になる x を求める：

sol2 = Solve[D[g, x] == 0, x]
これより sol2 として x = −1, 0 得る．そこでの g の値をg/.sol2として求めると，それぞれ g(−1) = 0, g(0) = −6 であることが分かる．最後に 2 階微分を計算させて，

4

15

章末「研究」のヒントと略解

D[g, x, x] /.sol2より g(−1) = −36, g(0) = 36．よって x = −1 で極大値 0，x = 0 で極小値

−6.
(4) 上の (1) と (3) を参照．
(5) r を変化させたとき，グラフを集合として{}
G(r) = (x, y, z) | z = x2 + y 2 , |x| ≤ r, |y| ≤ rと表してみる．この集合は縦・横・高さの比が 2r : 2r : 2r 2 の箱に入っているが，Mathematica はこれを縦・横・高さの比が 1 : 1 : 0.4 の箱に収まるように(X, Y, Z) = (ax, by, cz) の形の座標変換をして表示していると考えられる．
（Plot3D のオプション BoxRatios のデフォルト値は 1 : 1 : 0.4 だから．
）実際にそのような箱B = {(X, Y, Z) | |X| ≤ 1/2, |Y | ≤ 1/2, 0 ≤ Z ≤ 2/5}をとり，座標変換の係数を計算すると a = 1/(2r), b = 1/(2r), c = 1/(5r 2 ) となることがわかる．したがって関係式 z = x2 + y 2 は Z = 4(X 2 + Y 2 )/5 と同値になり，これはもとのグラフ G(r) に依存しない(6) f[x_,y_] = x^2 + y^2; r = 2;
Plot3D[ f[x,y], {x, -r, r}, {y, -r, r},
BoxRatios -> {1, 1, 1}, PlotRange -> {0, 4},
ClippingStyle -> None]

6.6 研究(1) 略．
(2) 土星：sph[s_, t_] = {Cos[s] Sin[t], Sin[s] Sin[t], Cos[t]};ParametricPlot3D[{sph[u, v], {(1.3 + v/5) Cos[u], (1.3 + v/5) Sin[u], 0}},{v, 0, Pi}, {u, 0, 2 Pi},
Mesh -> None, Boxed -> False, Axes -> False,
PlotStyle->{Opacity[0.9, Yellow], Opacity[0.7, Yellow]}](3) メビウスの帯：ParametricPlot3D[5{Cos[t] (3 + r Cos[t/2]), Sin[t](3 + r Cos[t/2]),r Sin[t/2] }, {r, -1, 1}, {t, 0, 2 Pi}]7.6 研究（その他の研究課題）(1) 好 み も あ る だ ろ う が ，た と え ば 問 題 7.2 で 定 義 し た 関 数 に 加 え て
pt[p_] := Graphics[{PointSize[Large], Red, Point[{p, f[p]}]}]
gr[p_] := Plot[{f[x], tanline[x, p]},{x, -3, 3},
PlotRange -> {-10, 10}, PlotStyle -> {Thick, Thin}];のように定義しておき，

Manipulate[ Show[{gr[p], pt[p]}], {{p, -1.5}, -3, 3}]としてみよう．次のような結果が得られる：

p105-3-21-123-5-10(2) z = f (x, y) = xy の接平面を図示せよ，という問題．多変数微分積分の講義でプレゼンする場面をイメージしている．まずは簡単に，次のようにすれば「動く接平面グラフ」が得られるだろう：

L = 5;
f[x_, y_] = x y;fx[x_, y_] = D[f[x, y], x];fy[x_, y_] = D[f[x, y], y];tanPl[p_, q_, x_, y_] = fx[p, q](x - p) + fy[p, q](y - q) + f[p, q];
ここまでが描画範囲を決める定数 L と接平面を表す関数の定義である．次に fのグラフを定義：

gr = Plot3D[f[x, y], {x, -L, L}, {y, -L, L},615

章末「研究」のヒントと略解

PlotStyle-> Directive[Opacity[0.6], Blue] ];接 平 面 が 隠 れ な い よ う に 透 明 度 を 高 め ，色 を 便 宜 的 に 青 に し た ．次 に

Manipulate で接平面を動かすための関数を定義：
pl[p_, q_] := Plot3D[tanPl[p, q, x, y],{x, p - L, p + L},{y, q - L, q + L}];
最後に Manipulate で動かす：
Manipulate[ Show[gr, pl[p, q]],{{p, 0}, -L, L},{{q, 0}, -L, L}]

それなりのものは得られるが，あまり見やすいグラフではない，と思われたのではないだろうか．
次に，見映えの良さを高めるために，（ちょっと大変かもしれないが）改めて以下のように入力してみよう．*2 まずは関数の定義：

L = 5;
f[x_, y_] = x y;fx[x_, y_] = D[f[x, y], x];fy[x_, y_] = D[f[x, y], y];tanx[p_, q_, t_] = fx[p, q](t - p) + f[p, q];tany[p_, q_, t_] = fy[p, q](t - q) + f[p, q];*2 見やすさには人それぞれ好みもあるだろう．あくまで参考までに．

7

tanPl[p_, q_, x_, y_] = fx[p, q](x - p) + fy[p, q](y - q) + f[p, q];
接線も一緒に描画するための関数も追加した．（偏微分係数の意味を伝えるため．）次に f のグラフ：

gr = Plot3D[f[x, y], {x, -L, L}, {y, -L, L},
PlotPoints -> 30,
PlotRange -> 2 L,
BoxRatios -> {1, 1, 1},
PlotStyle -> Directive[Opacity[0.8],
Specularity[White,50]],
Mesh -> None,
ClippingStyle -> None];

次に点 (p, q) における 2 本の接線と，それらの「影」を箱の底面に描くための関数を定義：

lin[p_,q_] := ParametricPlot3D[{{t, q, tanx[p, q, t]},{p, t, tany[p, q, t]},{t, q, -2 L}, {p, t, -2 L}},{t, -L, L}, PlotPoints-> 10,
PlotStyle -> {{Thickness[0.007], Blue},{Thickness[0.007], Red},Blue, Red} ];次に点 (p, q) における接平面とその「影」を箱の底面に描くための関数を定義：

pla[p_, q_] := Plot3D[{tanPl[p, q, x, y], -2 L},{x, p - L/2, p + L/2},{y, q - L/2, q + L/2},
PlotPoints -> 10,
Mesh -> None,
PlotStyle ->{Directive[Opacity[0.6], Green],Directive[Opacity[0.2], Green]}];
接平面だけは描画範囲を小さめに描画し，箱の底面にも「影」が描かれるようにした．あとは Manipulate で動かす：

Manipulate[ Show[gr, lin[p, q], pla[p, q],815

章末「研究」のヒントと略解

AxesLabel -> {x, y, z}],{{p, 0}, -L, L}, {{q, 0}, -L, L} ]以上のように入力すると，次のような出力が得られる：(3) 高次元のテイラー展開・フーリエ展開についても原理的には 7.4 節のような「動くグラフ」が描けるのだが，計算に時間がかかるので Manipulate の結果はスムーズに動作しないようだ．

8.7 研究(1) 略．
(2) まずは関数と変数（極座標）を定義
f[z_] := z^2; z := r Exp[I t];
たとえば次のように入力すると単位円板の f による像が得られる：

ParametricPlot[vect[f[z]], {r, 0, 1}, {t, 0, 2 Pi}]
円板の中心と半径を Manipulate で動かすための関数を定義*3 ：*3 8 章の [32] もこのように書き直したほうがわかりやすいかもしれない．9

mapf[w_, a_] := ParametricPlot[{vect[z+ w], vect[f[z+w]]},{r, 0, a}, {t, 0, 2 Pi}, PlotRange -> 5]
あとはManipulate[mapf[comp[p], a],{{p, {-2, -1}}, Locator}, {{a, 1.5}, 0, 3}]と入力すれば 110 ページの図のような結果が得られる．
(3) 正方形とその像を横に並べて表示する方法．f は上と同じだが，z の定義を変える：

r = 0.7; f[z_] := z^2; z := x + y I;
sqr[w_] := ParametricPlot[vect[w + z],{x, -r, r}, {y, -r, r}, PlotRange -> 5]
fsqr[w_]:= ParametricPlot[vect[f[w + z]],{x, -r, r}, {y, -r, r}, PlotRange -> 5]
Manipulate[{{sqr[comp[p]], fsqr[comp[p]]}}//TableForm,{{p, {-2, -1}}, Locator}]
正方形とその像を表示する関数を別に定義し，Manipulate でまとめて表示する．その際「TableForm を 1 × 2 行列に適用する」形を取った．（本のなかではこうした tricky な方法は避けたが，使える手法だとは思う．）
結果は次のようになる：

442200-2-2-4-4-4-2024-4-2024(4) たとえば次のように入力する:
L = 1.3 Pi; f[z_] := Exp[z];lines[{a_, b_}] := ParametricPlot[{vect[a + I t], vect[ t + I b]},{t, -L, L}, PlotRange -> L]flines[{a_, b_}] := ParametricPlot[

1015

章末「研究」のヒントと略解{vect[f[a + I t]], vect[f[t + I b]]},{t, -L, L}, PlotRange -> L]
Manipulate[{{lines[p], flines[p]}}//TableForm,{{p, {1, 1}}, Locator}]
こうするとつぎのような出力が得られる：

-44

4

2

22-24-42-2-2-2-4-44

ポ イ ン ト を 挙 げ る な ら ，直 行 す る 2 直 線 を 描 か せ る 関 数 lines をlines[a_, b_] ではなくロケータにあわせて lines[{a_, b_}] と定義して い る こ と（ そ の 像 flines も 同 様 ），(3) の よ う に Manipulate の 中 でTableForm を用いたことであろう．

11.10 研究√√(1) 訂正：本では bn = ϕn+1 / 5 (n ≥ 0) とありますが，ここでは bn = ϕn / 5 (n ≥1) としておきます．（得られる数列自体は同じ．）√√
方程式 x2 = x + 1 の解を ϕ = (1 + 5)/2, ψ = (1 − 5)/2 とおくと（ϕ は黄金比），フィボナッチ数列はFn =ϕn − ψ nϕn − ψ n√=ϕ−ψ5√と書けることがしられている．いま |Fn − bn | = |ψ n |/ 5 であるが，|ψ| < 1より Fn と bn の誤差は指数関数的に小さくなるのである．

Mathematica で確認するには，Table と TableForm で誤差を表にするのがよいだろう．ListLinePlot を使っても，ふたつの数列は近すぎて違いが表現できない．
（誤差を対数グラフで表現するのは面白いかもしれない．）(2) a[x_] := If[x < 0, 0, Exp[-1/x^2]]

11b[x_] := a[x]/(a[x] + a[1 - x])c[x_, e_] := b[x/e + 2]*b[-x/e + 2]
aa = Plot[a[x], {x, -1, 2}, AspectRatio->Automatic];
bb = Plot[b[x], {x, -1, 2}, AspectRatio->Automatic];
cc = Plot[c[x,0.5], {x, -2, 2}, AspectRatio->Automatic];{{aa,bb,cc}}//TableForm(3) 素数の組 (p, q, r)（ただし p < q < r ）が三つ子素数であれば，r = p + 6 である．よって次のようにすればよい：

condiTriple[{p_, q_, r_}] := (r - p == 6);
tab = Table[{Prime[n], Prime[n + 1], Prime[n + 2]},{n, 1, 1000}];
triples = Select[tab, condiTriple]
1002 番目の素数までに 99 組の三つ子素数が見つかるはずである．
(4) li = {PrimePi[x], x/(Log[x] - 1.08366), RiemannR[x]};{{Plot[li, {x, 0, 100}], Plot[li,{x, 0, 10000}]}}//TableFormと入力する．結果は次のようになる：
120040100030800206001040020020406080100-102000

40006000

8000 10 000

1 万以下ではルジャンドルもリーマンも非常に良い近似となっている．ガウス，チェビシェフとの比較も同様にやってみるとよい．

12.7 研究• 3 状態にしたければ，mod 2 のかわりに mod 3 とすればよい．例えば
g[x_, y_, z_] := Mod[x + y + z + y*z, 3];
色をつけるには，本文にある方法よりも手っ取り早く，ArrayPlot のオプショ
ン ColorFunction を使うのもよいだろう．たとえばArrayPlot[..., ColorFunction -> Hue]
ArrayPlot[..., ColorFunction -> "TemperatureMap"]のように関数名やカラースキーム名で指定すればよい．

1215• n = 101; g[x_, y_, z_] :=章末「研究」のヒントと略解

Mod[x + y + z + y*z, 2];としよう．まずは固定的境界条件．両端を 1 とした：

rule1[c_, i_]:= If[i == 1, g[1, c[[1]], c[[2]]],
If[i == n, g[c[[n-1]], c[[n]], 1],g[c[[i-1]], c[[i]], c[[i+1]]]]]
反射的境界条件：

rule2[c_, i_]:= If[i == 1, g[c[[2]], c[[1]], c[[2]]],
If[i == n, g[c[[n-1]], c[[n]], c[[n-1]]],g[c[[i-1]], c[[i]], c[[i+1]]]]]と入力する．結果は次のようになる：

左端が周期的，中央が固定的（両端は 1），右側が反射的．