Mathematica入門

Mathematica 9.0
2014/4/27 Takuichi Hirano (Tokyo Institute of Technology)

【はじめに】
・Mathematicaは記号数式処理が得意であり、それが他の言語との違いである。
・基本的にMathematicaはインタープリタ言語である
・実行はShift+Return
・Ctrl+kで補完
・乗算はスペース" "かアスタリスク"*"
・冪乗は"^"。ただし、数式入力パレットを使うと見やすい。除算やルート"Sqrt"なども同じ
・論理演算のandは"&&"、orは"||"。
・比較演算子の等号（イコール）は"=="。（=は代入演算子）
・組み込み関数は頭が大文字、それ以外は小文字になっている。?で組み込みかどうか調べることができる。組み込み関数との判別のために、自分が使う変数は頭を小文字にしておくとよい。変数名に"_"は使えない
・%は１つ前の実行結果。%%は２つ前の結果。%nはOut[n]の実行結果
・Mathematicaファイル(.nb)をノートブックと言う
・[]は関数定義用、{}はリスト用、()は演算順位指定用。
・「パレット->BasicMathInput」で数式入力ボタンを表示
・「評価->評価を中断」で評価を中断
・「評価->カーネルを終了->Local（ローカル）」でカーネル（バックグランドで動いている計算プロセス）を終了
・「評価->ノートブックを評価」でファイル全体を実行
・セル(右の"]")を選択し、「書式」でいろいろな書式設定
・「セル->全出力セルの削除」で出力セルを削除（保存するとファイルサイズが小さくなる）
・「(*」と「*)」で囲まれた部分はコメント

数値計算

電卓のように算術計算をすることもできる。

5+7

電卓とは違って、 Mathematica は完全に正確な結果を与える。

(3^100)

桁を保証する
N[]は有限の桁精度で表示するコマンド

N[Pi,100]

%    (* 1つ前の実行結果 *)

%%   (* 2つ前の実行結果 *)

%1    (* Out[1]の実行結果 *)

円周率

N[Pi]

N[π]    (* パレットにある *)

自然対数の底

N[E]

N[E]    (* パレットにある *)

変数

代入演算子(=)は右の値を左の変数に代入する。
改行によって複数のコマンドをまとめて実行できる。;もコマンドの区切りであり、;があると結果は画面に出力されない。

a=1;
b=3;

a/b

a=1.;

a/b

Clear[a,b]

a/b

関数の定義

x_の_は全てのパターンを意味する。
:=(Set Delayed)は=(Set)とは違って、右辺をすぐには評価しないで、使われるときになって初めて評価を開始する。
普通の関数定義は以下の表現を用いる。

f[x_]:=(1/(Sqrt[2 π]))*Exp[-((x^2)/2)]  (* 正規分布 *)

?f  (* 関数定義、変数の値などを表示 *)

g[x]の定義はg[x]にしかマッチしない。

g[x]:=x+1

g[2]

g[x]

g[x]の使用例。再帰呼出（数学用語で言う漸化式）の定義の終了条件など。

fact[0]=1;
fact[n_]:=n*fact[n-1];

fact[3]

グラフィック

関数 f[x] のグラフを描く。

f[r_]:=(1/r)

Plot[Abs[f[r]],{r,-3,3}]

グラフのフォーマット。
３次元グラフでもほぼ同じオプションが通用する。詳しくはヘルプで"Plot", "Plot3D"などで検索し、例を見て真似する。

Options[Plot]

Plot[Abs[f[r]],{r,-3,3},
PlotStyle->{Red},
Axes->None,
Frame->True,
PlotRange->{{-3,3},{0,10}},
FrameLabel->{"r","f(r)","",""}]

Plot[{1,x,(x^2),(x^3)},{x,-2,2}]

場合分けのグラフ

g[x_]:=x /; x≤0;
g[x_]:=0/; (x≥0) && (x≤1);
g[x_]:=(((x-1))^2) /; x≥1;

Plot[g[x],{x,-1,2}]

2 変数関数の 3D グラフを描く。

Plot3D[f[(Sqrt[(x^2)+(y^2)])],{x,-3,3},{y,-3,3}]

等高線

g1=ContourPlot[f[(Sqrt[(x^2)+(y^2)])],{x,-3,3},{y,-3,3},
ContourShading->None]

ベクトル図

<<****`はパッケージ****を読み込む命令。
"VectorFieldPlots"を「ヘルプ->ドキュメントセンター」で検索してもよいし、「ヘルプ->ドキュメントセンター->データの可視化」の「その他」の項目を見てもいろいろなパッケージがあるのがわかる。

<<VectorFieldPlots`

g2=VectorFieldPlot[(1/(Sqrt[(x^2)+(y^2)]))*(({x,y})/(Sqrt[(x^2)+(y^2)])),{x,-3,3},{y,-3,3}]

Show[{g1,g2},Graphics[{{RGBColor[1,0,1],Disk[{0,0},0.1]},Circle[{0,0},0.1]}]]

代数、微積分

代数

(a+b x+c y) (d+e x+f y)

展開する。

Expand[%]

因数分解する。

Factor[%]

極限を求める。

Limit[((Sin[x])/x),x->0]

Limit[((Cos[x])^(1/(x^2))),x->0]

Limit[(((1+(1/n)))^n),n->∞]
(* Definition of the base of the natural log *)

三角関数

Sin[a]+Sin[b]

TrigFactor[%]

TrigExpand[%]

TrigReduce[%]

TrigToExp[%]

Simplify[%]

FullSimplify[%]

級数 (テイラー, ローラン展開)

0 を中心とした第 5 項目までのテイラー展開(マクローリン展開)。

Series[Exp[x],{x,0,5}]

1 を中心とした第 4 項目までのテイラー展開。

Series[Log[x],{x,1,4}]

I を中心とした第 5 項目までのローラン展開。

Series[1/(1+z^2),{z,I,5}]

微分

(∂/(∂x))(Sin[x] Cos[2 x] Exp[3 x]+(x^3))

D[Sin[x] Cos[2 x] Exp[3 x]+(x^3),x]

微分方程式

DSolve[y'[x]==a y[x],y[x],x]

ベッセルの微分方程式
特殊関数もいろいろある

DSolve[(x^2)*D[y[x],{x,2}]+x*D[y[x],x]+((x^2)-(α^2))*y[x]==0,y[x],x]

解の利用方法
"expression /. {x->y,z->w}"は左の表現expressionにおいて表現xをyに、zをwに変換する。

sol=DSolve[y''[x]+(k^2)*y[x]==0,y[x],x]

f[x_]=(y[x] /. sol[[1]] /. {C[1]->0,C[2]->1})

Plot[{f[x] /.k->1,f[x] /.k->2,f[x] /.k->3},{x,0,2*π}]

積分

∫(Sqrt[1+(x^2)]) +(x^2) Sin[x] Cos[2 x] dx

Integrate[Sqrt[1+(x^2)]+(x^2) Sin[x] Cos[2 x],x]

Integrate[x*(y^2)*Sin[x]*Cos[x],{x,a,b},{y,c,d}]

∫((Sin[x])/x) dx
不定積分はできない（初等関数およびその四則演算と合成関数では表現できない）が、特殊関数が定義されている。

Integrate[((Sin[x])/x),x]

数値積分

Integrate[x*(y^2)*Sin[x]*Cos[x],{x,π,2*π},{y,1,2*π}]

NIntegrate[x*(y^2)*Sin[x]*Cos[x],{x,π,2*π},{y,1,2*π}]

線形代数

ベクトル、行列
MathematicaではListとして扱われる。

ベクトル（Mathematicaの内部処理用の用語は「リスト」。数学のベクトルに限定しない。）

v1={x1,y1,z1};
v2={x2,y2,z2};

v1[[1]]    (* 成分の抽出 *)

v1.v2    (* 内積 *)

Cross[v1,v2]    (* 外積 *)

Drop[v1,-1]    (* 最後の要素の削除 *)

Drop[v1,{2}]    (* ２番目の要素の削除 *)

Take[v1,{2,3}]    (* ２番目から３番目までの要素を抽出 *)

Append[v1,w1]    (* 要素を追加 *)

Join[v1,v2]    (* リストを結合 *)

行列（２次元リスト）

m1={{a1,b1},{c1,d1}};
m2={{a2,b2},{c2,d2}};

m1[[1,2]]    (* 要素の抽出 *)

MatrixForm[m1]    (* 行列形式で表示 *)

MatrixForm[*] と * // MatrixForm は同じ意味。後から処理を加えたいときに便利な記法。 * // Simplify などもよく用いる。

m2 // MatrixForm    (* 行列形式で表示 *)

m1+m2 // MatrixForm    (* 行列の加算 *)

m1*m2 // MatrixForm    (* 行列の要素同士の乗算 *)

m1.m2 // MatrixForm    (* 行列の積 *)

m1.Drop[v1,-1] // MatrixForm    (* 行列とベクトルの積 *)

Transpose[m1] // MatrixForm   (* 転置 *)

規則あるリストの生成

Table[(i^2),{i,1,5}]

Table[(i^2)*j,{j,1,3},{i,1,5}]

リストへの関数の適用

Clear[f,lis];

lis=Table[(i^2),{i,1,5}]

Map[f,lis]

Apply[f,lis]

行列 m とベクトル v の定義。
BasicMathInputのパレットの(({{□, □}, {□, □}}))を使うこともできる。Ctrl+Returnで行を追加、Ctrl+,で列を追加。

m=(({{a, b}, {c, d}}));
v=(({{x}, {y}}));

行列とベクトルの掛け算。

m.v

m.Inverse[m]

簡単な表現に直す。

Simplify[%]

行列式を求める。

Det[m]

逆行列を求める。

Inverse[m]

特性多項式を求める。

CharacteristicPolynomial[m,x]

方程式を解く。

Solve[%==0,x]

固有値を求める。

Eigenvalues[m]

固有ベクトルを求める。

Eigenvectors[m]

方程式を解く

bについて解く

Solve[a*x+b*y==c,b]

解の公式

Solve[a*(x^2)+b*x+c==0,x]

Solve[a*(x^3)+b*(x^2)+c*x+d==0,x] // Simplify

Solve[a*(x^4)+b*(x^3)+c*(x^2)+d*x+e==0,x]

5次以上の方程式には解の公式が存在しない（ガロア理論）。つまり、解を初等関数とその合成関数および四則演算で表すことができない。もちろん、解は複素数の範囲には存在する（代数学の基本定理）。

Solve[a*(x^5)+b*(x^4)+c*(x^3)+d*(x^2)+e*x+f==0,x]

Solve[{x+2 y==1,5 x+4 y==1}]

数値的に解く

NSolve[{x+2 y==1,5 x+4 y==1},{x,y}]

Solveコマンドは非線形連立方程式も扱うこうとができる

Solve[{x+2 y==1,5 x+4 (y^2)==1},{x,y}]

ベクトル解析

カルテシアン(x,y,z)座標系

Grad[f[x,y,z],{x,y,z},"Cartesian"]

Grad[x*(y^2)*(z^3),{x,y,z},"Cartesian"]

Div[{fx[x,y,z],fy[x,y,z],fz[x,y,z]},{x,y,z},"Cartesian"]

Curl[{fx[x,y,z],fy[x,y,z],fz[x,y,z]},{x,y,z},"Cartesian"]

ベクトル公式の確認

Div[Curl[{fx[x,y,z],fy[x,y,z],fz[x,y,z]},{x,y,z},"Cartesian"],{x,y,z},"Cartesian"]

Curl[Grad[f[x,y,z],{x,y,z},"Cartesian"],{x,y,z},"Cartesian"]

Div[Grad[f[x,y,z],{x,y,z},"Cartesian"],{x,y,z},"Cartesian"]

-Curl[Curl[{fx[x,y,z],fy[x,y,z],fz[x,y,z]},{x,y,z},"Cartesian"],{x,y,z},"Cartesian"]+Grad[Div[{fx[x,y,z],fy[x,y,z],fz[x,y,z]},{x,y,z},"Cartesian"],{x,y,z},"Cartesian"]    (* Vector Laplacian *)

円筒(ρ,φ,z)座標系

Grad[f[ρ,φ,z],{ρ,φ,z},"Cylindrical"]

Div[{fρ[ρ,φ,z],fφ[ρ,φ,z],fz[ρ,φ,z]},{ρ,φ,z},"Cylindrical"]

Curl[{fρ[ρ,φ,z],fφ[ρ,φ,z],fz[ρ,φ,z]},{ρ,φ,z},"Cylindrical"]

球(r,θ,φ)座標系

Grad[f[r,θ,φ],{r,θ,φ},"Spherical"]

Div[{fr[r,θ,φ],fθ[r,θ,φ],fφ[r,θ,φ]},{r,θ,φ},"Spherical"]

Curl[{fr[r,θ,φ],fθ[r,θ,φ],fφ[r,θ,φ]},{r,θ,φ},"Spherical"]

複素関数論

複素数

z=1+I*2

Re[z]

Im[z]

Abs[z]

Conjugate[z]

ComplexExpand[Sin[x+I*y]]    (* 変数を全て実数と見なし、a+I*bの形にする *)

ComplexExpand[x*Sin[y],{x,y}]    (* x,yは複素数とし、他の変数を全て実数と見なし、a+I*bの形にする *)

他にも"ComplexMap"をヘルプで検索すると複素平面を描く面白い例を見ることができる

留数

Clear[z]

Res[(1/(z^2))+(2/z)+3+4 z,0]: 極の例 (ローラン展開の1/zの項の係数に等しい)

Residue[(1/(z^2))+(2/z)+3+4 z,{z,0}]

Res[((Sin[z])/z),0]: 除去可能な特異点の例 (一見特異点のようだが、実際には有限値なので、1/zの項の係数は0である)

Residue[((Sin[z])/z),{z,0}]

Res[(e^(1/z)),0]: 真性特異点の例 (ローラン展開でも表現できない)

Residue[Exp[1/z],{z,0}]

フーリエ変換

記号的にフーリエ変換を行う
FT: F(ω)=(1/(Sqrt[2 π]))(∫_(-∞)^∞)f(t) (e^(-j ω t)) dt
IFT: f(t)=(1/(Sqrt[2 π]))(∫_(-∞)^∞)F(ω) (e^(j ω t)) dω

FourierTransform[1,t,ω]

FourierTransform[UnitStep[t+τ/2]-UnitStep[t-τ/2],t,ω]

DFT: X(k)=(1/(Sqrt[n]))*(∑_(j=1)^n)x(j)*Exp[-2*π*(k-1)*(j-1)/n]
IDFT:  x(j)=(1/(Sqrt[n]))*(∑_(k=1)^n)x(j)*Exp[-2*π*(k-1)*(j-1)/n]
ω=2π (k-1)/(N Δt), t=(j-1) Δt

num1=10;
lis=Join[Table[1,{i,1,num1}],Table[0,{i,1,40}]]
n=Length[lis];
Δt=0.5;
flis=Δt*(Sqrt[n])*Fourier[lis];
ListPlot[lis]
g1=ListLinePlot[Abs[flis],PlotRange->All]
g2=ListLinePlot[Table[{2*π*(k-1)/(n*Δt),Abs[flis][[k]]},{k,1,n}],PlotRange->All]

τ=num1*Δt;
g3=Plot[Abs[((Sqrt[2*π]))*(((Sqrt[(2/π)]) Sin[((τ ω)/2)])/ω)],{ω,0,2*π/Δt},PlotRange->All,PlotStyle->Red]
Show[g2,g3]

手続き型プログラミング

繰り返し(Do)

If文とBreakでループの途中中断可。
If文とContinueでループ内の残りの行を飛ばし、ループの次のステップに進む。

Moduleで手続きを関数化できる。

和を求めるプログラムの例

s=0;
Do[
s+=i;
If[Mod[i,10]==0,Print["i=",i]];
,{i,1,100}
]
Print["Sum=",s];