1 今回使うプログラムの基礎
コンピュータプログラムでは
• 条件文
• ループ文
の2つが使えればなんとかなるものです．今回はそれらについて述べましょう．
1.1 条件文
Ifを用いる文です．条件をみたしているか，いないかによって，以降で違 う働きをします．
Mathematicaでは，最も簡単な形は If[条件, 真の場合の操作]
です．例えば，
If[x>0, y=1]
ならば，xが正のとき，yを1にします．同じようなものでは If[x>0, y=1,y=-1]
とすると，xが正のとき，yを1に，xが負のときには，yを−1にします．つ まり
If[条件, 真の場合の操作，偽の場合の操作]
では，xが正のとき，yを1に，xが負のときには，yを−1に，さらにxが 0ならばyも0にするにはどうしたらいいでしょう．真が偽かの判定は1つ しかできませんから，2つIf文を使わなければなりません．
If[x>0, y=1,If[x==0,y=0,y=-1]]
とすればよいのです．ここで，気をつけなければいけないのは，等しいのは
“==”とすることです．コンピュータのプログラムではA=BとはBの内 容をAにしまうという代入文に使われるので，等しいを表すには==としま す．この他に<=とか>=も使えます．意味はわかりますね．この他に!=と いうのがあって，等しくないという意味です．
条件をみたしたときに，複数の操作をさせましょう．例えば，xが正なら ば，yを1にzを2にするには，真の場合の動作2つを1つの文で表せばよ いだけですから，それは
y=1;z=2
と間をセミコロンでつなぐと1つの文とみなせますので If[x>0,y=1;z=2]
とすればよいのです．
それでは複数の条件に対応するときにはどうでしょう．それには
• 条件Aと条件Bを同時にみたしているか，A&&B
• 条件Aか条件Bをみたしているか，A||B の2種類があります．
xもyも正のときに，zを1にするには If[(x>0)&&(y>0), z=1]
xかyが正のときに，zを1にするには If[(x>0)||(y>0), z=1]
Mathematicaには便利な関数があって，
• 偶数であるかどうかを判断するEvenQ
• 奇数であるかどうかを判断するOddQ
• 素数であるかどうかを判断するPrimeQ
• 整数であるかどうかを判断するIntegerQ
などたくさんあります．これらは，真ならTrue，偽ならFalseを返してきま す．したがって，xが偶数ならばyは1，奇数ならば−1にするには
If[EvenQ[x],y=1,y=-1]
とすればよいのです．
問題1 ガウス記号を表すFloor関数を用いて，EvenQとIntegerQを作って ください．
解
If[Floor[x/2]*2 == x, Print[True], Print[False]]
If[x - Floor[x] > 0, Print[False], Print[True]]
1.2 ループ文
繰り返しをさせる文で，もっともコンピュータらしい文です．いくつかの タイプがありますが，ここではFor文を学びましょう．
For[初期条件,終了条件,継続文，実行文]
の形をしています．継続文とは実行を終わったときに次に何をするかを決め る文です．1から10まで数字を印刷させるには
For[i=1,i<=10,i++,Print[i] ]
ここで，継続文は i++
ですが，これはiに1を加えるという意味です．同様なやり方として i+=1
とか i=i+1
としても同じです．コンピュータ特有の表現ですが，
For[i=1,i<=10,i*=2,Print[i] ]
の意味はわかりますね．1から10まで足すプログラムならば，足し算をした 結果をしまう場所をsとして，初めはs=0として
s=0;For[i=1,i<=10,i++,s+=i ]
これだけでは，結果がでないので，終わったらsを表示するように s=0;For[i=1,i<=10,i++,s+=i ];s
とすればよいわけです．また，sの初期値を設定するのをFor文の中にして For[i=1;s=0,i<=10,i++,s+=i ];s
偶数だけを足すにはどうしたらよいでしょう．最も簡単なのはiの初期値を2 にして，2ずつ加えて
For[i=2;s=0,i<=10,i+=2,s+=i ];s
もしくは足すときに偶数であるかどうかを判断して，
For[i=1;s=0,i<=10,i++,If[EvenQ[i],s+=i] ];s 問題2 100までの素数の積を求めてください．．
解
For[i=1;s=1,i<=100,i++,If[PrimeQ[i],s*=i]];s
問題3 PrimeQを作ってください．
解
For[i = 2; s = 0, i < x, i++, If[Floor[x/i]*i == x, s++]];
If[ s == 0, Print[True], Print[False]]
ここでsは約数の個数になっていますし，
For[i = 2; s = 0, i < x, i++, If[Floor[x/i]*i == x, s++;Print[i] ] ] とすると約数が表示されます．
1.3 ソートプログラム
データx1, x2, . . . , xnを大きさの順に並べ直すことをソートすると言いま す．このプログラムを作ってみましょう．
データを名前をつけてListにします．例えば L=List[2,3,1,7,9,8]
リストの2番目はL[[2]]で表せます．まず，1番大きいのを1番から6番ま で探します．それには，1番目とi番目を比べて，i番目の方が大きかったら，
入れ替えればよいので
For[i=2,i<=Length[L],i++,If[L[[i]]<L[[1]],L[[1]]=L[[i]] ] ];L
ここでLength[L]はリストLの長さ，上の例なら6になります．これを実行 すると
L={1, 3, 1, 7, 9, 8}
と確かに1番目は最小値になりますが，1番目のデータ2が消えてしまいま す．それでは入れ替えはできません．i番目が大きかったら，1番目とi番目 を入れ替える必要があります．
For[i=2,i<=Length[L],i++,If[L[[i]]<L[[1]],L[[1]]=L[[i]];L[[i]]=L[[1]] ] ];L
でいいでしょうか．やってみればわかりますが，入れ替わっていません．ど うしてでしょうか．最初にL[[1]]=L[[i]]を実行していますから，既に1番に はi番目のデータが入っています．その次にL[[i]]=L[[1]]とやっても，1番に はもともとはi番のデータが入っているのですから，何も行われないことに なってしまうのです．それを解消するには1番目のデータを別に保存をして おいて，
For[i=2,i<=Length[L],i++,If[L[[i]]<L[[1]], tmp=L[[1]];L[[1]]=L[[i]];L[[i]]=tmp ] ];L
とすればよいことになります．2番目に小さなデータは2番目と3番以降を 比較すればよいので
For[i=3,i<=Length[L],i++,If[L[[i]]<L[[2]], tmp=L[[2]];L[[2]]=L[[i]];L[[i]]=tmp ] ];L
と順番にLength[L]−1番までやれば，よいのでこれもループ文にすると
For[j=1,j<Length[L],j++,
For[i=j+1,i<=Length[L],i++,If[L[[i]]<L[[j]], tmp=L[[j]];L[[j]]=L[[i]];L[[i]]=tmp ] ] ];L
これで確かに並べ替えはできますが，データ数が多くなるととても時間がか かって実用にはなりません．様々な工夫が行われていて，バブルソート，シェ ルソート，クイックソートなどが有名です．クイックソートは，軸になる数 を決めて，それより大きいグループと小さいグループに分けることを繰り返 す方法です．実際に100枚程度の並べ替えでも上の方法では大変ですが，ク イックソートをするとかなり速くなります．
上のデータなら
1. 最初の軸を3.5にすると，{2,1,3}と {7,9,8}
2. 最初のグループの軸を2.5，後ろのグループの軸を8.5とすると，{2,1}, {3}, {7,8},{9}となります．
3. 最後に1番目のグループの軸を 1.5，３番目のグループの軸を8.5とし て順番に並びます．
2 Newton 法
一般的に方程式f(x) = 0を解くには，
Solve[f[x]==0,x]
のようにすれば解けます．しかし，ご存知のように5次方程式以上では一般 的に解く方法はありません．数値解を求めるには
NSolve[f[x]==0,x]
とやると，多項式の解ぐらいなら求まりますが，exなどが入ると実用にはな りません．
そこで，一般的な方程式f(x) = 0の解を漸近的に求める方法です．収束が 速く実用になります．実はMathematicaではFindRootとして実装されてい ます．
FindRoot[f[x]==0,{x,x0}]
とします．x0は計算を始める値です．複数の解があるときには初期値ごとに 解が異なります．
Newton法は
1. f(x0)における接線を求める．
2. 接線がx軸と交わる点をx1とおく．
3. このx1を初期値として，1を実行する．
この繰り返しです．どこで止めるかを決めておかないといつまでも計算し続 けます．
x0での接線は
y=f0(x0)(x−x0) +f(x0) です．したがって，x軸との交点は
x1=−f(x0) f0(x0)+x0 まず，f を入力し，その微分を
Df[x_]=D[f[x],x]
で求めておきます．初期値もxとしておきます．n回で計算を打ち切るならば f[x_]=E^x-x-2;
Df[x_]=D[f[x],x];
x=1;
n=10;
For[i=1,i<=n,i++, x=- f[x]/Df[x]+x ];x
とやればいいのですが，これだと式をそのまま実行するので遅くなります，そ こで，小数計算にするように
f[x_]=E^x-x-2;
Df[x_]=D[f[x],x];
x=1;
n=10;
For[i=1,i<=n,i++, x=N[-f[x]/Df[x]+x]
];x
とやれば，x= 1.14619になります．実は既に3回目でこの値になっています．
3 Taylor 展開
関数fのaにおけるTaylor展開とは，収束半径をρとすると，|x−a|< ρ ならば
f(x) =f(a) +f0(a)(x−a) +f00(a)
2 (x−a)2+f000(a)
3! (x−a)3+· · · です．Mathematicaには
Series[Sin[x],{x,1,3}]
のような関数があって，sinxをa= 1中心のテイラー展開の3次までを求め てくれます．このプログラムを作って，身近な関数を近似してみましょう．
ex = 1 +x+x2 2 +x3
3! +· · · cosx = 1−x2
2 +x4 4! −x6
6! +· · · sinx = x−x3
3! +x5 5! −x7
7! +· · · log(1−x) = −x−x2
2 −x3 3 − · · · arctanx = x−x3
3 +x5 5 −x7
7 +· · · 1
1−x = 1 +x+x2+· · · µ 1
1−x
¶2
= 1 + 2x+ 3x2+ 4x3+· · ·
√1 +x = 1 +1 2x+1
2· −1 2 ·x2
2 +1 2 ·−1
2 ·−3 2 ·x3
3! +· · · 問題4
arctan 1 = 1−1 3 +1
5 +· · ·
を用いて，πが3.14であることをチェックするにはいくつまで計算すればよ いでしょうか．
解 294まで計算すると四捨五入して，3.14になります． いくらでも加え ていくのでは，コンピュータが止まりません．いつ止めたらいいでしょう．そ れをするためには，誤差項が必要なのです．n−1次までで近似したとき f(x) =f(a)+f0(a)(x−a)+f00(a)
2 (x−a)2+· · ·+f(n−1)(a)
(n−1)! (x−a)n−1+Rn(x) と表せて， Rn(x)が誤差項です．誤差項にはいろいろな表現がありますが，
もっとも一般的なのはaとxの間のある点cが存在して Rn(x) =f(n)(c)
n! (x−a)n です，
例えば，sin 1を誤差が 10001 で求めたいならば，f(n)は ±sinか ±cosな ので，絶対値は1以下ですからa= 0として
|Rn(x)|<|x|n n!
を得ます．x= 1のときを求めるのですから，6! = 720,7! = 5040ですから，
n= 7とすればよいことになります．こうして止めるnが定まれば，後はプ ログラムです．
まず，関数fのk回微分は D[f[x],{x,k}]
です．これのx=aでの値は，xに aを代入するには D[f[x],{x,k}]/.{x->a}
とします．これをk!で割るには，代入が終わってから，k!で割るので，括弧 でかこって
(D[f[x],{x,k}]/.{x->a})/k!
とします．
f[x_]=Sin[x];n=7;
For[k=1;s=f[0],k<=n,k++,
s+= (D[f[x],{x,k}]/.{x->0})/k!]
;s
とすればよいのですが，これだと，式のまま計算してしまうので，小数計算 をするようにするには，初期値を小数にしておけばよいので
f[x_] = Sin[x]; n = 8;
For[k = 1; s = N[f[0]], k <= n, k++, s += (D[f[x], {x, k}] /. {x -> 0})/k!]; s となります．一般にa中心のテイラー展開にするには f[x_] = Sin[x]; n = 8;a=0;
For[k = 1; s = N[f[0]], k <= n, k++, s += (D[f[x], {x, k}] /. {x -> a})/k!]; s
さらに，x= 1でなくても使えるようにするには，上のxをyに取り替えて f[y_] = Sin[y]; n = 8;a=0;x=1;
For[k = 1; s = N[f[0]], k <= n, k++,
s += (D[f[y], {y, k}] /. {y -> a})*x^k/k!]; s これだと，sin 2も上の式でx= 2とすれば出せますが，
|Rn(2)| ≤ 2n n!
ですから ，29!9 = 0.0014,210!10 = 0.00028ですので，この場合には n= 10にし なければならないのです．
グラフで比較するならば，x= 1を除いて f[y_] = Sin[y]; n = 8;a=0;
For[k = 1; s = N[f[0]], k <= n, k++,
s += (D[f[y], {y, k}] /. {y -> a})*x^k/k!]; s としてから
Show[Plot[s, {x, 0, 3Pi/2}], Plot[f[x], {x, 0, 3Pi/2}, PlotStyle -> {Hue[0]}]]
とすると，sinxの収束半径は∞ですから，とてもよく近似していることが わかります．しかし，例えば√
1 +xのように収束半径が!だと，−1< x <1 の範囲で収束はしますが，±1を越えて近似はできません．
問題5 cos 3， π，eを誤差 10001 で求めるにはどうしたらよいか，プログラ ムを作ってください．
解 誤差項をきちんと求めること，eの計算では Rn(x) =ec
n!
を考慮して，e <3を用いるとよい．
このnも定めるようにプログラムを改良しましょう．sinxのa= 0における テイラー展開のn次の誤差は |x|n!n より小であることから
f[y_] = Sin[y]; a=0;x=1;gosa=1/1000;
For[n=1,Abs[x]^n/n!>gosa,n++,];
For[k = 1; s = N[f[0]], k <= n, k++,
s += (D[f[y], {y, k}] /. {y -> a})*x^k/k!]; s
誤差がgosa以下になるようなnを定めるだけですので，実行文はありません．
4 Digital 技術の基礎， Fourier 級数
計量線形空間V の正規直交基底e1,e2, . . .とは，長さが1で直交すること です．内積を用いると
(ei,ei) = 1, (ei,ej) = 0 (i6=j) をみたすことです．周期2πをもつ関数の空間では
Z 2π
0
f(x)g(x)dx= (f, g)
を内積とすることができます．ここでg(x)はg(x)の複素共役です．より正 確に述べるならば，|f(x)|2が積分可能な複素関数全体は上の内積により計量 線形空間(L2[0,2π]という)になります．ただし，次元は可算無限次元です．
この空間は量子力学でも用いられる重要な空間です．内積からはノルム(つ まり，ベクトルの長さ)が
||f||=p
(f, f) = Z 2π
0
|f(x)|2dx
で定められます．ここで，f(x)f(x) =|f(x)|2であることに注意しておきま す．このノルムから，関数f とgとの距離が
d(f, g) =||f −g||
によって定められます．この距離でL2[0,2π]は位相空間となります．
この基底(正規ではない)として
1,cosx,cos 2x, . . . ,sinx,sin 2x, . . .
を選ぶことができます．これらが直交することは，nとmにさまざまな値を 入れて
Integrate[Cos[n*x],{x,0,2Pi}]
Integrate[Sin[n*x],{x,0,2Pi}]
Integrate[Cos[n*x]*Cos[m*x],{x,0,2Pi}]
Integrate[SIn[n*x]*Sin[m*x],{x,0,2Pi}]
Integrate[Cos[n*x]*Sin[m*x],{x,0,2Pi}]
を実行することにより確かめられます． 正規直交基底にするには長さを1に すればよいので
Sqrt[Integrate[1^2,{x,0,2Pi}]]
Sqrt[Integrate[Cos[n*x]^2,{x,0,2Pi}]]
Sqrt[Integrate[Sin[n*x]^2,{x,0,2Pi}]]
を行うことで
√1 2π, 1
√πcosx, 1
√πcos 2x, . . . , 1
√πsinx, 1
√πsin 2x, . . . が正規直交基底になることがわかります．
したがって，任意の関数f ∈L2[0,2π]はこれらの和で f(x) =a0 1
√2π+a1 1
√πcosx+a2 1
√πcos 2x+· · ·+b1 1
√πsinx+b2 1
√πsin 2x+· · · と表されることになります．ただし，収束するというのはn次まで加えた fn(x) =a0 1
√2π+a1 1
√πcosx+· · ·+an 1
√πcosnx+b1 1
√πsinx+· · ·+bn 1
√πsinnx がf(x)に近づくことではなく，fnと f との距離||fn−f||が0に近づくこ とです．
関数fに対して，a0, a1, . . . , b1, b2, . . .をfのフーリエ級数と言います．逆 に数列a0, a1, . . . , b1, b2, . . .が与えられれば，上の式で関数fを求めることが できます．実際の計算では，双方に√
πをつけると煩わしいだけでなく，計 算誤差も増えるので，関数から級数を求めるときか，級数から関数を求める ときのどちらかにπをつけます．これがデジタル技術の基礎理論です．
数列から関数を求める場合，任意の 数列でよいわけではありません．f ∈ L2[0,2π]ですから，fの長さ||f||は有限でなければなりませんので，
||f||2= (f, f) = Z 2π
0
|f(x)|2dx= X∞
n=0
|an|2+ X∞
n=1
|bn|2
ですから，右辺が有限なもの(l2に属するという)でなければなりません．こ のことは，量子力学が作られた時代に，無限次元のベクトルとして，その上 の無限次元の行列を用いて構築する一派(ハイゼンベルグ)と関数の上に偏微 分方程式を用いて構築する一派(シュレディンガー )が互いの主張をしあっ て，相手を非難したのですが，実は同じものだったという，今では笑い話が あるのです．
これらのanやbnを決めるには，例えば (f, 1
√πcosnx) = (a0 1
√2π+a1 1
√πcosx+a2 1
√πcos 2x+· · · +b1 1
√πsinx+b2 1
√πsin 2x+· · · , 1
√πcosx)
= an1
π(cosnx,cosnx) =an
などとすればよいわけです．
まず，
s = Integrate[f[x]/Sqrt[2 Pi], {x, 0, 2 Pi}]/Sqrt[2 Pi]
とすれば，定数に対応する項が求まります．残りはリストにしましょう．n 項まで求めるならば
a = Table[ Integrate[f[x]*Cos[k*x]/Sqrt[Pi], {x, 0, 2 Pi}], {k, 1, n}];
b = Table[Integrate[f[x]*Sin[k*x]/Sqrt[Pi], {x, 0, 2 Pi}], {k, 1, n}];
とすればよいことになります．これらを加えればよいのですから For[k = 1, k <= n, k++,
s += a[[k]]*Cos[k*x]/Sqrt[Pi] + b[[k]]*Sin[k*x]/Sqrt[Pi]];
となります．最後に，元の式と比較する図を描きます．全体をまとめると n = 10;
f[x_] = x;
s = Integrate[f[x]/Sqrt[2 Pi], {x, 0, 2 Pi}]/Sqrt[2 Pi];
a = Table[
Integrate[f[x]*Cos[k*x]/Sqrt[Pi], {x, 0, 2 Pi}], {k, 1, n}];
b = Table[
Integrate[f[x]*Sin[k*x]/Sqrt[Pi], {x, 0, 2 Pi}], {k, 1, n}];
For[k = 1, k <= n, k++,
s += a[[k]]*Cos[k*x]/Sqrt[Pi] + b[[k]]*Sin[k*x]/Sqrt[Pi]];
Show[Plot[s, {x, 0, 2 Pi}],
Plot[f[x], {x, 0, 2 Pi}, PlotStyle -> {Hue[0]}]]
n= 10程度でもかなり近いことがわかりますが，端の方は近づいていないこ とがわかります．例えば周期 π2 の鋸波
f(x) =x (mod π 2) とか方形波
f(x) =
1 0≤x < π2, π≤x < 32π
−1 その他
を描いてみましょう．これらはテレビやコンピュータなどで重要な役割を果 たす関数ですので，これらを正確に作ることは電子工学において欠かせない ものです．これらを描いてみると，フーリエ級数での近似は，各点収束はし ていないことが見えてきます．まして，一様収束などはしません．収束の概 念ができていない頃はこのことがわからなくて，大変苦労をしたようで，こ のことがかえって，計算方法の進歩をもたらした歴史があります．
問題6 鋸波，方形波について実行してみてください．
解
f[x_] = Mod[x, Pi/2];
f[x_] = If[(x < Pi/2) || ((x > Pi) && (x < 3 Pi/2)), 1, -1];
