MATLAB Overview

• 「MATLABハンドブック」，小林一行，秀和
• 「MATLABグラフィクス集」，小国力，朝倉書店
• 「MATLABと利用の実際---現代の応用数学とCG」， 小国力 ，サイエンス社
• 「Introduction To Scilab の日本語訳」， http://www.geocities.jp/rui_hirokawa/scilab/ （Scilab 日本語ページ）
• 「Scilab 入門(2002年度版)」，高橋義男（沼津高専教授）， 沼津高専 電気電子工学科４年生の講義用手作資料また，MATLABのリンクは，目次のところにまとめました。
  

  ---

  実行例

  ここに，MATLABを実際に使って簡単な計算をした例を示します。
  図の後に，各部を説明します。
  

  ◎MATLAB Command Window [＿][□][×] File Edit Window Help [New M-File] [Open File] . . . >> a=4; b=5*a, c=b-10 - - - - - - - (1) b = 20 c = 10 >> x=[0 1 2 3]*30 - - - - - - - (2) x = 0 30 60 90 >> y=sin(x*pi/180) - - - - - - - (3) y = 0 0.5000 0.8660 1.0000 >> - - - - - - - (4) ◎Scilab [＿][□][×] File Edit Contral Demos GraphicWindow0 Help Editor --> a=4; b=5*a, c=b-10 - - - - - - - (1) b = 20. c = 10. --> x=[0 1 2 3]*30 - - - - - - - (2) x = ! 0. 30. 60. 90. ! --> y=sin(x*%pi/180) - - - - - - - (3) y = ! 0. 0.5 0.8660254 1. ! --> - - - - - - - (4)

  ここでは，MATLAB実行画面を示します。 本来はもっと空白行が多いですが，ここでは，空白行の一部は詰めました ここでは入力はこの色で示します。(本来はモノクロ画面) (括弧)はコメント用に私が書き加えたものです (1) : 対話式に入力します。 普通のコンピュータ言語(BASIC, C など)のように式を書けます。 なお，BASICと同様，定義しない変数を使えます。 直ぐに計算し，結果を表示します。 ";"や","で区切れば1行に複数の式も書けます。 (";"区切りは，結果表示無し) (2) : MATLABの特徴である，配列の一括処理の例です。 0 から 3 までの配列を定義し，要素それぞれに 30 を掛けます。 計算結果も，要素それぞれを示します。 (3) : pi (=π=3.14･･･) / 180 の掛け算も，sin という計算も， 配列全体に一括処理できます。 30度刻みで sin(0度)から，sin(90度)まで計算出来ました。 (4) : 入力プロンプトが出て，次の入力を待っています。 以下，同様に処理を続けられます。 Scilab の場合，見てくれが多少変わります。 行列の表記が明示的になりました。 数値の表示法も変わりました。 円周率(=pi) は MATLAB と同じく最初から組み込まれていますが， MATLAB では "pi" という変数だったのに対して Scilab では "%pi" です。 プロンプトは " --> " になりました。 メニューが変わりました。

  
  

  ---

  MATLAB環境の説明(ディレクトリ，ファイル，m-ファイル)

  MATLABでは，演算結果をファイルに保存することができます。
  現在使用中のディレクトリはカレントディレクトリです。
  m-ファイルとは，拡張子.m がつくテキスト形式のファイルで， MATLABのプログラムを入れておくことができます。
  (Scilab では，.sce や .sci ファイル になります)
  代表的には以下のようなコマンドを覚えておけば大丈夫。
  

  MATLAB環境 (Scilabに違いがあれば)

  helpdesk MATLAB 全般に関するヘルプ。メニューのヘルプからでも可能 （Scilab ではメニューの help から） lookfor keyword keyword に関連するコマンドを挙げる （Scilab では apropos keyword） help keyword keyword に関するマニュアル カーソルキー ↑で，前に入力した行の復帰。何行でも戻れるし，再編集できます。 demo デモンストレーション起動 （Scilabではメニューから起動） exit MATLABの終了 ls または dir カレントディレクトリのファイルのリストを見る （Scilabでは unix_w dir と打ち込む。 unix コマンドは，OS (含 Win.) のコマンドを実行する。 _w は，出力を Scilab ウインドウ内に設定） pwd カレントディレクトリのパスを見る cd path カレントディレクトリを，path に設定する （Scilab では chdir path） edit エディタが起動。これで新たな m-ファイルを編集できる。 edit test test.m という m-ファイルを編集 type test test.m という m-ファイルの中身を表示 delete test.m test.m という m-ファイルを消去 （Scilab では .sce ファイルを使います。 詳しくは「プログラミング」の節で）

  
  

  ---

  演算a --- 電卓として使う，関数一覧

  MATLABの最も簡単な応用方法は，電卓として利用することです。
  倍精度型が標準ですから，高い精度の計算が保証されます。
  変数を使えは，計算結果を一時的に記憶しておくことができます。 その場合，「変数 = 数式」 という書き方をするのは，他の言語と全く同じです。 この「変数 = 数式」という形式を代入文と呼びます。 MATLABでは，前もって定義すること無しに変数を使えます。 なお，一時的というのは， MATLABを終了したり，clear コマンドを実行したり，同じ変数に対して改めて代入するまでという意味です。
  代表的には以下のようなコマンドや，使い方を覚えておけば大丈夫。
  

  演算---電卓として使う，関数一覧 (Scilabに違いがあれば)

  演算と代入文 例:"y = y + 1" y+1 という演算をし，その結果を y に代入。表示もする 変数のみ + (C/R) 例:"a" 変数 a の内容を表示させる ";" 例:"y = a + b ;" 最後のセミコロンは，セパレータで，かつ，"表示無し"の意味 ans "a+b" のように演算だけを実行させると，デフォルト変数"ans"に答えが代入される 四則演算とべき乗 a + b, a - b, a * b, a / b (or b \ a), a ^ b "help ==" 左のように打てば，演算子のリストを見ることができる （Scilabの似たコマンドは調査中） inf 1/0 のこと:無限大 NaN または nan Not a Number : 数字じゃない pi 円周率 （Scilabでは %pi と表記） i または j 虚数 (i^2 = j^2 = -1) 例:"a=1-2i" （Scilabでは %i と表記） save myworkspace myworkspace.mat というファイルに，使用中の変数全体を保管する ( Scilabは拡張子を補わないので， この例の通りに打つと myworkspace というファイルに保管する） load myworkspace myworkspace.mat というファイルから，変数を復帰する （Scilabは拡張子を補わないので， この例の通りに打つと myworkspace というファイルから復帰する） who または whos 使用中の変数一覧を見る。whoは簡略版，whosは詳細版 （Scilab には who コマンドのみ存在し，それが詳細版） clear 変数をすべて破棄する clear a a という変数を破棄する 宣言 MATLAB は，C や Pascal と異なり，宣言しない変数を使うことができる。 変数名は，1文字目がアルファベットで， 2文字目以降は数字とアンダースコア'_'も許される。 A と a MATLAB は大文字と小文字を区別する。A と a は別物であり，同時に共存できる （Scilab の最新版も同じく大文字と小文字を区別する。なお， Scilabの古いバージョンでは大文字と小文字を区別しない） format 数値の表示形式を標準とする （Scilab では format('v',10) ） format long format long e format long g 固定小数点15桁表示 浮動小数点15桁表示（Scilab では format('e',15) ） 固定or浮動小数点15桁表示（Scilab では format('v',15) ） なお，内部の計算は常に有効桁数約16桁の浮動小数点形式で行われている。 formatコマンドで表示形式を変わっても内部の数値が変わるわけではない。 sin, sqrt など y=sin(x) など，普通の関数は普通にあります("help sin" など打ち込んで確かめよう) y = real(v * exp(j*2*pi*f*t)) 殆んどの関数も複素演算できるので，電気回路理論と整合が良い

  
  

  ---

  演算b --- 1次元配列，配列のグラフ化

  MATLABの良いところは，配列を簡単に扱える点です。
  配列を使って代入文を書くと，自動的に"要素ごとの演算"をやってくれます。
  以下にその利点を例を示します。

  MATLAB	BASIC (MATLABに対応した書き方)	BASIC (BASICに最適化)
  >> x=linspace(0,pi,101); >> y=sin(x); >> plot(x,y,'b');	dim x(101),y(101) for i=0 to 100 x(i)=pi*i/100 next i for i=0 to 100 y(i)=sin( x(i) ) next i ;グラフの座標系の定義 ;グラフの縦軸，横軸，目盛を書く for i=0 to 99 line(x(i),y(i))-(x(i+1),y(i+1)) next i	;グラフの座標系の定義 ;グラフの縦軸，横軸，目盛を書く for i=0 to 100 x=pi*i/100 y=sin( x ) if i>0 then line(x0,y0)-(x,y) x0=x y0=y next i

  上記の例では MATLAB ではずいぶんと簡潔に書けています。 しかも標準のサブルーチン plot のおかげで， 簡単にグラフが描けてしまいます。 ここで BASIC の名誉のために言っておきますが， MATLAB と BASIC のどちらも， もう一方に対して劣っているわけではありません。 両方の言語ともに適した分野が有り， たまたま上記のような例は， MATLAB が最も得意としたということです。 すなわち，配列を使い，その結果を効率良く表示するという分野です。
  ところで，MATLABでは，１次元配列を扱うとき， 関数の引数等でちょっと不思議なところがあります。（例 ： size 関数） これは，MATLAB がもともと MATRIX (2次元配列) を扱うことを前提としているからでしょう。
  代表的には以下のようなコマンドや，使い方を覚えておけば大丈夫。
  

  演算---1次元配列，配列のグラフ化 (Scilabに違いがあれば)

  x(3) x という配列の3番目の要素 x=[2 7 3] x(1)=2, x(2)=7, x(3)=3 という配列を作る x(2)=4 y=x 配列の要素に対して代入できる。 配列全体も式一つで代入できる x=[a b] 配列の連結 (5:1:12) [5 6 7 8 9 10 11 12] "5 から始まり，1 づつ増加し，12 までの配列" (5:12) (5:1:12)のこと (5:2:10) [5 7 9] "5 から始まり，2 づつ増加し，10 までの配列" a(5:2:10) [a(5) a(7) a(9)] linspace(0,pi,101) 等差数列：0 から pi までを，101個の要素で logspace(2,5,10) 等比数列：10^2 から 10^5 までを，10個の要素で + - マトリクス演算 == 要素ごとの演算 * / マトリクス演算 ~= 要素ごとの演算 .* ./ .\ .^ マトリクス演算 ~= 要素ごとの演算 2次元のグラフ figure(n) n という番号のついたグラフ用紙を用意します。 デフォルトで figure(1) が指定済み。 いったん figure(n) と宣言をすると， グラフ関連のコマンドは figure(n) に反映します。 （Scilab のグラフック窓関連は調査中） plot(x,y) 座標(x(n), y(n)) (n=1,2,...,配列の最後まで)をつないだグラフ plot という関数は，3つ以上のパラメタも指定できる。詳しくは， "help plot" と打ち込んで学ぼう。そのヘルプでは， loglog，semilog，極座標(polar) など別形式のグラフも紹介している （Scilab の plot 関数はパラメタが 2 つの場合は同様のグラフを書くが， 3 つ目以上のパラメタは MATLAB とは異なる。 なお，loglog のグラフについては， plot2d 関数のところで述べる。） 同じ座標に3つのグラフを表示 最初に次の操作により 3 つの関数を定義する。 これらを一つグラフに書いてみよう。 x1=0:0.1:7; y1=sin(1*x); y2=sin(2*x); y3=sin(3*x); MATLAB では，デフォルト設定では plot 関数が呼び出されると 先ず画面を消してからグラフを書く。 そのため，そのままでは 3 つのグラフを画面上に書くことは出来ない。 hold on という設定により画面を消さなくなり， hold off でデフォルトに戻る。 なお，plot は一口に解説しきれないくらいの機能を持っているので， help plot と打って勉強して欲しい。 上記の plot 関数の3番目の引数は線の色を指定するものである。 MATLAB plot(x,y1,'b') hold on plot(x,y2,'r') plot(x,y3,'g') hold off Scilab xbasc() plot2d(x,y1,style=1) plot2d(x,y2,style=2) plot2d(x,y3,style=3) （Scilab では，xbasc により画面を消去する。 plot は簡単な機能しか持たないようであるが， plot2d はそこそこの能力を持っているので plot2d を使う。 plot2d は前の画面を消さずにグラフを書く。 （だから xbasc が必要） plot2d の style パラメタは色を指定するが， 負の値を設定するとプロットになる。 なお，plot2d の拡張形は plot2d('nn',x,y,style=n) となる。 最初の 'nn' は両軸とも通常目盛りを意味する。 デフォルト設定でもあり指定無き場合は'nn'とみなされる。 代わりに 'ln' とすれば横軸だけが Log であり， 'nl' なら縦軸だけが Log であり， 'll' なら両対数グラフとなる。 大文字による指定（'NL'など）は禁止である） これ以降，配列専用の関数-1 (1次元でも良く使うもの) max, min, sort など max(a) とは a という配列の最大値を求める （1次元配列については，普通に最大値が計算されるが， 2次元配列に対して max を適用すると， 各行の最大値が計算されて，1次元配列が返ってくる） （Scilab では max(2 次元配列) などの結果はスカラー値になる） length(a) a が m行×n列 の 2次元配列であるとし， mとnのうちの大きな方の数が返ってくる。 （Scilab では 要素の数，即ち m×n が返ってくる） size(a) a という配列のサイズを求める。 ただし，2次元配列が仮定されているので， なぜか 何行×何列 という答えが返ってくる その他 上記の例はすべて 1行×n列 の配列。 要素をセミコロンで区切ることにより (例 a=[1;2;3;4;5])， n行×1列の配列も作ることができる

  
  操作の例：

  CD に刻まれる音を題材にしてグラフ表示をしてみよう。 ここでは，1 ms の間の現象を確認したいと考える。 MATLAB では取り扱うべき変数 t を，Δt という微小時間ごと取り扱う。 本来，時間は連続的に変化するが， ここでは Δt を十分に小さくすることで良しと考える。 >> t = 0:1E-7:0.001; >> size(t) 注意：1E-7 とは，1×10-7という意味である。即ち，0.1μs毎に扱う。 ここで，1kHz = 1000 Hz の正弦波を考える。 >> f1 = sin(2*pi*1000*t); これをプロットしてみる。 >> plot(t,f1); カーソルを使えば元の行に直ぐに戻れる。 2kHz，10kHz，50kHz など幾つかの関数をプロットしてみよう。 以上は特に CD に限らない話であった。 ところで CD には，44.1 kHz という有限の時間でサンプリングされた情報が刻まれている。 （ここではグラフの書きやすさからサンプリング周波数を 40 kHz として話を進める。） すなわち，時間は飛び飛びになるのである。 fclk = 40 kHz ということは， Tclk = 25 μs である。 先ほどの例で言えば，t を 250 個ごとに使うことに対応する。 それでは，サンプリングした情報をグラフ化してみよう。 >> plot(t(1:250:10001),f1(1:250:10001),'o'); 注意：plot 関数の 3 番目の引数の 'o' は，値をプロットすることを意味する 以上の操作のうち，MATLAB に打ち込む行を整理： >> t = 0:1E-7:0.001; >> size(t) >> f1 = sin(2*pi*1000*t); >> f1 = sin(2*pi*2000*t); >> f1 = sin(2*pi*10000*t); >> f1 = sin(2*pi*39000*t); >> f1 = sin(2*pi*41000*t); >> plot(t,f1); >> plot(t(1:250:10001),f1(1:250:10001),'o'); あるいは次の行を実行しておくと， >> tt = t(1:250:10001); あとは，次の 1 行で全てをまかなうことができる。 >> f=1000; plot(t,sin(2*pi*f*t),'r',tt,sin(2*pi*f*tt),'o'); 注意：MATLAB の plot 関数をうまく使うと，1 行で 2 つ以上の関数を表示できます。 Scilab 用に打ち込む行を整理： --> t = 0:1E-7:0.001; --> size(t) --> f1 = sin(2*%pi*1000*t); --> f1 = sin(2*%pi*2000*t); --> f1 = sin(2*%pi*10000*t); --> f1 = sin(2*%pi*39000*t); --> f1 = sin(2*%pi*41000*t); --> plot2d(t,f1); --> plot2d(t(1:250:10001),f1(1:250:10001),-1); あるいは次の行を実行しておくと， --> tt = t(1:250:10001); あとは，次の 1 行で全てをまかなうことができる。 --> f=1000; xbasc(); plot2d(t,sin(2*%pi*f*t)); plot2d(tt,sin(2*%pi*f*tt),-1);

  
  

  ---

  演算c --- 文字列

  文字列演算は，コンソール出力を行う際に必要な処理です。
  代表的には以下のようなコマンドや，使い方を覚えておけば大丈夫。
  

  文字列 ---きれいな出力を得るために--- (Scilabに違いがあれば)

  MATLABの場合 概要 MATLABでは，文字列は 1行×n列 の配列として取り扱います。 代入 例:" msg='television' " 文字列はシングルクォートで囲む 連結 例:" b=['This is a ' msg '!!!'] " 配列の連結と同じである 抜き出し 例:" c=msg(3:5) " 3文字目から5文字目まで(上の例からでは 'lev') 数値<==文字 例:" a = 5 + str2num('123') " 文字<==数値 例:" s = ['The answer is ' num2str(a * b)] " disp(a) disp関数の引数は一つしか許されない。 disp(a) では a という変数の内容を表示。 （単に a(c/r) では 'a=' という表示もついてくる） a は，文字列・数値・配列のいずれでも可能 応用例 例:" disp( ['Data = ' num2str(x) ', Its squear =' num2str(x * x)] ) " Scilab の場合 (文字列の扱い法は MATLAB とは大きく異なる) 概要 Scilab では，文字列は専用の文字列方のようです。 代入 Scilab でも，" msg='television' " 連結 Scilab では， b = 'This is a ' + msg + '!!!' 抜き出し Scilab を調査中 数値<==文字 Scilab を調査中 文字<==数値 Scilab では， s = 'The answer is ' + string(a * b) 応用例 Scilab例:" disp( 'Data = ' + string(x) + ', Its squear =' string(x * x) ) "

  
  

  ---

  演算d --- 2次元配列，配列のグラフ化

  当然，2次元配列を扱うのも MATLAB の得意とするところです。
  要素の取り扱い方法など，MATLAB 固有の簡潔な表記法があります。
  代表的には以下のようなコマンドや，使い方を覚えておけば大丈夫。
  

  演算---2次元配列，配列のグラフ化 (Scilabに違いがあれば)

  概要 1次元配列の場合の演算はたいていそのまま使える ";"を使った入力 例:" x = [1 2 3;4 5 6]" とは， x = 1 2 3 4 5 6 要素 a(1,2) とは 2 のこと : a の1行目でかつ2列目の要素 a(2,:) とは [4 5 6] のこと : a の2行目 a(:,2) とは [2;5] のこと : a の2列目 a(:,1:2) とは [1 2;4 5] のこと a(:,1:2:3) とは [1 3;4 6] のこと 行と列の入替 a=x' 連結 x=[a b] : 列方向(右方向)に連結 x=[a;b] : 行方向(下方向)に連結 変数を 2 つ持つ 3 次元グラフ （MATLABの場合） 3Dグラフ [x,y]=meshgrid(-2:0.5:2, -2:0.1:2); z=cos( x.^2 + y.^2 ); mesh(x,y,z); この場合，meshgrid により x と y に変数域を設定する。 x も y も 9列(x軸)41行(y軸) のマトリクスになる。 x は x 座標(横軸)に対してのみ変化。 z=cos(...) は演算の例。z も xやyと同じサイズのマトリクス。 mesh 関数は 3 次元グラフを描く。 "help mesh" を見ると，その他の3次元プロットも紹介してくれる 視点 view([-1 -3 5]); どの方向からグラフを見るか視点を指定できる 変数を 2 つ持つ 3 次元グラフ （Scilabの場合） 3Dグラフ （Scilab の場合も plot3d(x,y,z) によってメッシュ表示をしてくれる。 z が二次元配列（マトリクス）である点は MATLAB と同様であるが， Scilab の場合は x と y として一次元配列（ベクトル） を使う点が MATLAB と異なる。 以下のように演算すれば，MATLAB の例と同様のグラフを描ける。 x=-2:0.5:2; y=-2:0.1:2; xmat=x'*ones(1,length(y)); ymat=ones(length(x),1)*y; z=cos( xmat.^2+ymat.^2 ); xbasc(); plot3d(x,y,z); なお，plot3d も plot2d 同様に，そのままでは画面を消去しないため， ここでは xbasc() ... 画面消去 ... を加えた。 視点 Scilab を調査中 これ以降，配列専用の関数-2 (2次元配列) zeros zeros(4,4) とは 4行4列 の，zeros(3,5) とは 3行5列の，全ての要素が 零 の配列 （zeros(4,4) のように 2 つの引数が一致するなら， MATLAB では zeros(4) という形式に略しても同じだが， Scilab では zeros(4) ではマトリクスは作れない） ones zeros と同様だが，全ての要素が 1 の配列 rand zeros と同様だが，それぞれの要素が 0 から 1 の間の乱数から成る配列 length(a) 配列 a について，行と列の長さの長いほうを与える （Scilab では 要素の数，即ち 行の数×列の数 が返ってくる）

  
  

  ---

  演算e --- マトリクス

  MATLAB の名前から予期できる通り，マトリクス演算は MATLAB が最も得意とするところです。 （もちろん Scilab も得意としています。）
  代表的には以下のようなコマンドや，使い方を覚えておけば大丈夫。
  

  演算---マトリクス (Scilabに違いがあれば)
  概要 配列の演算はたいてい使える 行列の設定 zeros(3,3) や，ones(3,3) や rand(3,3) などは， それぞれ 零，1，乱数 の 3行3列の行列である。 この表記から判るように，マトリクスを扱うのは得意である 演算 a=[1 2 3;4 5 6 7 8 0] b=[366; 804; 351] det(a) で，a の行列式が求められる。 x=inv(a)*b または x=a\b で，行列式" a x = b "を満たす x が求められる その他 a' とは，a の行と列を入れ替えたもの

  
  

  ---

  プログラミング(m-ファイル)

  一連の処理をテキストファイルにして拡張子に .m を付けておくと， そのファイル名を打つことによって一連の処理を自動的に実行してくれます。 拡張子が .m のファイルを m-ファイルと呼びます。
  （Scilab では .sce ファイルです。ただし .sci も使えるようです）
  代表的には以下のようなコマンドや，使い方を覚えておけば大丈夫。
  

  プログラミング(m-ファイル) (Scilabに違いがあれば)

  プログラム実行の概要（MATLABのプログラム環境） pwd ワーキングディレクトリの Path を表示 dir ワーキングディレクトリ内のファイルを表示 edit エディタが起動。これで新たな m-ファイルを編集できる。 edit mytest ワーキングディレクトリ内の mytest.m を編集 mytest ワーキングディレクトリ内の mytest.m というファイルに書かれたプログラムを実行させる （Scilab によるプログラム実行の概要 ... 拡張子は .sce） pwd ワーキングディレクトリの Path を表示 unix_w('dir') ワーキングディレクトリの中身を表示 エディタの起動 ツールバーの右端の Editor をクリックすると起動 なお，ver.2.7 よりも前の Scilab にはエディタはありません。 そういう時は，別のエディタでプログラムを書いてきました。 例えば，Windows に標準装備のメモ帳（notepad）でも十分実用に足ります。 フリーのエディタ（例えば Terapad ）を使うことで， 更に便利に編集できます。 ファイルの編集 エディタのメニューバーから編集するファイルを選ぶ exec mytest.sce ワーキングディレクトリ内の mytest.sce というファイルに書かれたプログラムを実行させる m-ファイルの作り方（MATLABプログラムの書き方） 概要 普通の演算はそのまま書いておけば大丈夫である コメント 行のうち % という文字以降はコメントとみなされる （Scilab では，// という２文字を使用） 一時停止 pause; キーを押すまで一時停止 （Scilab では， -1-> というプロンプトが現れて停止する。 継続には return を， 終了には abort を打ち込む。 pause(3); 3秒停止 （Scilab では，引数の形式は許されない） 入出力 a=b+c; 値をあえて表示する必要が無い時 a=b+c デバッグのためにいちいち値を知りたいとき disp(a) 変数名の表示は不用の時 disp(['a = ' num2str(a)]) 細かい表示付き （Scilab では文字列処理方法が異なる） a=input('Please input : a =') 数値入力 a=input('Please input : a =','s') 文字列入力 For for n=1:10 x(n)=sin(n*pi/10); end forのほかの書き方 for n=[4 5 9 10] for n=4:-1:-10 While while n<10 x(n)=sin(n*pi/10); n=n+1; end If if n==2 a=1; end if n==2 a=1; elseif n==3 a=9; else a=0; end Break break 文により，for または while ループから 抜け出すことができる。通常次のように使われる if n==2; break; end 関数 下のように，ファイル名と同じ名前の定義行が 1行目にある m-ファイルを用意しておくと， その関数が使えるようになる。 たとえば，下記の myfunc1.m が ワーキングディレクトリ内に存在すると，例えば x = myfunc1(4) によりこの関数が呼び出され，x は 24 になる。 ＜myfunc1.m＞ (Download) function y=myfunc1(x) s=1; for n=2:x s=s*n; end y=s; （Scilab による関数の取り扱いの概要 ... 拡張子は .sce） 関数 まず下のようなファイルをワーキングディレクトリ内に作る。 なお，関数の名前（プログラムの 1 行目参照）と， ファイル名は別で構わない。 続いて， getf('myfunc2.sce'); というコマンドを実行すると，ファイル myfunc2.sce が読み込まれ，関数 myfunc1 が使えるようになる。 例えば y = myfunc1(5) によりこの関数が呼び出され， y は 120 になる。 ＜myfunc2.sce＞ (Download) function y=myfunc1(x) s=1; for n=2:x s=s*n; end y=s; endfunction

  
  

  ---

  応用1a --- 交流回路の電圧波形(計算とグラフ化)

  MATLAB は複素計算もグラフ表示も簡単にできるため， 交流回路への応用も簡単です。
  実行例(他にも良いやり方はいろいろあります)
  

  応用---交流回路の電圧波形を計算(グラフの書き方) (Scilabに違いがあれば)

  系: R-L 直列回路 抵抗 コイル 入力電圧 定数値 200 Ohm 1 H 141 V 「50 Hzでの各部の波形をグラフ化するプログラム」 % prog1a1 - waveform of an R-L circuit R=200; L=1; V=141; f=50; t=linspace(0,0.05,1000); ZL=i*2*pi*f*L; Vin=V*exp(i*2*pi*f*t); VR=R/(R+ZL) * V*exp(i*2*pi*f*t); VL=ZL/(R+ZL) * V*exp(i*2*pi*f*t); plot(t,real(Vin),'g',t,real(VR),'r',t,real(VL),'b') % end of program // prog1a1s (Scilab) R=200; L=1; V=141; f=50; t=linspace(0,0.05,1000); ZL=%i*2*%pi*f*L; Vin=V*exp(%i*2*%pi*f*t); VR=R/(R+ZL) * V*exp(%i*2*%pi*f*t); VL=ZL/(R+ZL) * V*exp(%i*2*%pi*f*t); xbasc(); plot2d(t,real(Vin),style=1) plot2d(t,real(VR),style=2) plot2d(t,real(VL),style=3) 「各部の電圧の周波数依存性を求めるプログラム」 % prog1a2 - frequency responce of the circuit R=200; L=1; V=1; f=logspace(0,4,100); ZL=i*2*pi*f*L; VR=V * R ./ (R+ZL) ; VL=V * ZL ./ (R+ZL) ; loglog(f,abs(VR),'r',f,abs(VL),'b') % end of program // prog1a2s (Scilab) R=200; L=1; V=1; f=logspace(0,4,100); ZL=%i*2*%pi*f*L; VR=V * R ./ (R+ZL) ; VL=V * ZL ./ (R+ZL) ; xbasc(); plot2d('ll',f,abs(VR),style=1) （Scilab では plot2d で loglog のグラフも書ける。 plot2d の引数の最初の 'll'（小文字の L が二つ）は， 縦軸も横軸も Log であることを表す。 'ln' は横軸だけが Log であり，'nl' は縦軸だけが Log である。 'nn' はデフォルト設定であり，両軸とも通常目盛りを意味する。 本来は，VR と VL の2つのグラフを同時に書かせたいが， グラフの座標が変わってしまうので，ここでは一つずつ書かせる。 → VL も書かせてみよう % prog1a3 - frequency responce R1=1.e3; C1=1.e-6; V=1; f=logspace(-1,5,300); ZC=1 ./ (i*2*pi*f*C1); VR=V .* R1 ./ (R1+ZC) ; VC=V .* ZC ./ (R1+ZC) ; loglog(f,abs(VR),'r',f,abs(VC),'b') % end of program // prog1a3s (Scilab) R1=1.e3; C1=1.e-6; V=1; f=logspace(-1,5,300); ZC=1 ./ (%i*2*%pi*f*C1); VR=V .* R1 ./ (R1+ZC) ; VC=V .* ZC ./ (R1+ZC) ; xbasc(); plot2d('ll',f,abs(VR)) 注意：ある数 (例えば f ) の逆数を求める際は，空白を入れながら書こう。 (例えば 1 ./ f のように。) 空白無しに (1./f) と書いたならば， (1. / f) = (1.0 / f) = マトリクスの割り算 と判断されるのか，又は， (1 ./ f) = 要素ごとの割り算 と判断されるか分らない。 注意： V はスカラーであるため（配列ではないので）， 続く掛け算は * でも良いが， .* でもエラーにならないので， むしろ打ち間違いを防ぐために積極的に .* を使うことを勧める。 % prog1a4 - frequency responce R1=1.e3; R2=100.e3; C1=0.1e-6; C2=1.e-6; V=1; f=logspace(-1,5,300); ZC1=1 ./ (i*2*pi*f*C1); ZC2=1 ./ (i*2*pi*f*C2); Z2=R2 + ZC2; ZP=ZC1 .* Z2 ./ (ZC1 + Z2) ; VR1=V .* R1 ./ (R1 + ZP) ; VC1=V .* ZP ./ (R1 + ZP) ; VR2=VC1 .* R2 ./ Z2 ; VC2=VC1 .* ZC2 ./ Z2 ; loglog(f,abs(VR1),'r') hold on loglog(f,abs(VR2),'b') loglog(f,abs(VC1),'g') loglog(f,abs(VC2),'k') hold off % end of program // prog1a4s (Scilab) R1=1.e3; R2=100.e3; C1=0.1e-6; C2=1.e-6; V=1; f=logspace(-1,5,300); ZC1=1 ./ (%i*2*%pi*f*C1); ZC2=1 ./ (%i*2*%pi*f*C2); Z2=R2 + ZC2; ZP=ZC1 .* Z2 ./ (ZC1 + Z2) ; VR1=V .* R1 ./ (R1 + ZP) ; VC1=V .* ZP ./ (R1 + ZP) ; VR2=VC1 .* R2 ./ Z2 ; VC2=VC1 .* ZC2 ./ Z2 ; xbasc(); plot2d('ll',f,abs(VR)) 注意： 並列は，a. * b. / (a + b) である。

  
  

  ---

  応用1b --- 実験データの整理(データと，回帰線)

  MATLAB を使えば，電卓感覚で実験のデータ整理ができます。
  実行例(他にも良いやり方はいろいろあります)
  

  応用---実験データの整理(データと，回帰線) (Scilabに違いがあれば)

  例として，次の実験データについて回帰線を求める。 x(cm) 0 1 2 4 5 y(cm) 0.1 0.9 2.8 9.0 15.2 % prog1b1 - 回帰線 % 数値の設定 x=[0 1 2 4 5] y=[0.1 0.9 2.8 9.0 15.2] % polyfit により最少2乗法で1次式を求めます。 % 求まった係数を変数 p に代入します。 % もし数字 1 の代わりに 2 や 3 を入れれば， % その次数の回帰線が求まります [p,q]=polyfit(x,y,1) % グラフの横軸を設定します。 xx=0:0.1:5; % xxの各要素に対するyyの各要素を計算します。 yy=polyval(p,xx) % データは丸で，回帰線は青線を指定。 plot(x,y,'o',xx,yy,'g'); % end of program // prog1b1s (Scilab) x=[0 1 2 4 5] y=[0.1 0.9 2.8 9.0 15.2] [p,q]=reglin(x,y) xx=0:0.1:5; yy=xx*p + q; xbasc(); plot2d(x,y,style=-1); plot2d(xx,yy,style=2); （prog1b1s で reglin(x,y) は1次式近似を求めます。 polyfit で 1 を指定した場合と同じです。 Plot2d では style が負のときはプロットです。 ）

  
  

  ---

  応用1c --- Waveファイルと信号処理(FFT，フィルタ)

  MATLAB は信号処理も簡単に行なえます。
  実行例(他にも良いやり方はいろいろあります)
  

  応用---Waveファイルと信号処理(FFT，フィルタ) (Scilabは調査中)
  % prog1c1 - FFT % y にはサンプルした音の各時間ごとの強度が入ります。 % Fs には1秒あたりのサンプル数が， % Bits には信号のビット数が入ります。 % また信号がモノラルかステレオかにより % y は1行n列か，または2行n列になります。 [y,Fs,Bits]=wavread('sound1.wav'); % y のうち，2秒目から 4096 点を抜き出します。 % もしもCDと同じ44.1kSpsだったならば， % 4096 点とは 約0.09秒です yy=y(1,2*Fs:2*Fs+4095); % yy を再生し，表示します。 playsnd(yy,Fs,Bits); figure(1); plot((2*Fs:s*Fs+4096)/Fs, yy); % yy を sound2.wav というファイルに保存します。 wavwrite(yy,Fs,Bits,'sound2.wav'); % zz は yy のFFT結果が入ります。 % yy が 4096 = 2^12 点 なので，高速にフーリエ変換出来ます。 % zz もやはり 4096点 の配列になります。 % サンプリング周期 Fs と，4096点であることから， % abs(zz(2048))というデータは，周波数 Fs/2 [Hz] の % 周波数成分の強度を表わします。 zz=fft(yy); % zz は複素数ですから，そのままではグラフ化出来ないので， % 絶対値を取ってから表示します。 figure(2); plot((1:4096)*Fs/4096, abs(zz)); % end of program
  % prog1c2 - Filter % y には音データが入ります。 [y,Fs,Bits]=wavread('sound1.wav'); % 1000 Hz から 5000 Hz に通過域を持つ % 6 次のバタワース型フィルタを設計します。 % ナイキスト周波数を単位として周波数を設定します。 % "help butter"とすれば，他のフィルタの情報も得られます。 [b,a] = butter(6,[1000/(Fs/2) 5000/(Fs/2)]); % 前行で得たパラメタ a, b を使って % y をフィルタにかけ，出力 z を得ます。 z = filter(b,a,y); % end of program

  
  

  ---

  応用1d --- 数値積分

  数値積分も MATLAB は色々な関数を取り揃えてあり，簡単に行なえます。
  数値積分は，微分方程式の最も特殊な場合と考えることも出来ます。
  実行例(他にも良いやり方はいろいろあります)
  

  応用---数値積分 （Scilabに違いがあれば）

  まずは，以下の関数を作り，myfunc.m という名前で保存（Save)しておく。 これでファンクションが定義できたことになり，以降たとえば， y=myfunc(3.14159/2) とでも打ち込むと，y にはほぼ 1 (= sin(90度) )が設定される。 実際のところ積分の計算には必ずしもfunctionが必須というわけではないが， こうすることによって(1)プログラムが見やすくなる (2)function を使うと大きな利点もある(最後の例を参照のこと)ので，お勧めである（というか，functionの利用を前提として説明する）。 ＜myfunc3.m＞ (Download) function y=myfunc3(x) y=sin(x); % end of function （Scilab） ＜myfunc3.sce＞ (Download) function y=myfunc3(x) y=sin(x); endfunction さて，まずは自分でプログラムを作り， 区間 0〜pi にわたって積分しよう。 採用するアルゴリズムは，一番簡単な台形法とする。 T=h*(y1 + 2*y2 + 2*y3 + yn)/2 % prog1d1 - 台形法 n=201; x=linspace(0,pi,n); h2=(x(2)-x(1))/2; y=myfunc3(x); s=h2*(y(1)+2*sum(y(2:n-1))+y(n)) % end of program // prog1d1s (Scilab) // clear myfunc3 getf('myfunc3.sce'); n=201; x=linspace(0,%pi,n); h2=(x(2)-x(1))/2; y=myfunc3(x); s=h2*(y(1)+2*sum(y(2:n-1))+y(n)) （prog1d1s の clear 文は，メモリ内に残る myfunc3 を消す。 これは，myfunc3 を別の関数に変えたときに， メモリ内に残る myfunc3 を消すためである。 なお，clear 文の前の // は最初から消してしまってかまわない。 というのも， myfunc3 の読み込み前に clear myfunc3 を行ってもエラーにならないからである。） 上記の台形法と同じアルゴリズムで計算する関数が用意されており， 次のように行っても同じ結果が出る。 % prog1d2 - 積分の関数(アルゴリズムは台形法)を利用 n=201; x=linspace(0,pi,n); y=myfunc3(x); s=trapz(x,y) % end of program （Scilab については調査中） functionの形で積分関数が設定されている場合， MATLAB の組み込み関数を使うと，積分はたったの一命令で実行できる。 なお，'lookfor integral'で確認すると， 他の同様の関数を知ることが出来る。 組み込み関数による積分 % prog1d3 - integral s=quad('myfunc3',0,pi) % end of program // prog1d3s (Scilab) s=intc(0,%pi,myfunc3) 注意: MATLAB と Scilab を比較すると， 積分を行う関数は，名前が違うだけでなく， パラメタを与える順番も，パラメタの与え方 （'...クォーテーションマーク で囲むかどうか）まで異なる。

  
  

  ---

  応用1e --- 微分方程式

  微分方程式を考える上で，2次の微分方程式は基本中の基本で， 応用範囲の広いものです。 それは，変位と速度から加速度を導く系になぞらえることが出来ます。 また，RLC回路の応答にもなぞらえることが出来ます。
  そういうわけで，ここでは，2次の微分方程式の解法を例に取り， MATLAB で解く方法を考えたいと思います。
  取り扱う系は，振り子の振動を表わす，
  x'' = - a x' - sin(x)
  とします。
  実行例(他にも良いやり方はいろいろあります)
  

  応用---微分方程式 （Scilabに違いがあれば）

  ここでは，次の式を MATLAB を利用して解き，振り子の振動を解析するものとする。 x'' = - a x' - sin(x) 解説 ： もともとの系は振り子である。 変数 x [m] は振り子が動く経路上の長さを表す。 記号 ' は時間による微分を表す。 質量 m [kg] の質点が， 長さ L [m] の棒の先にくっついており， 空気抵抗は a [kg/s] ， 重力によって下の点に引き戻されようとする力は 重力加速度 g [kg m/s2] とsinを掛けた値になるというものであり， m x'' = - a x' - g sin(x/L) で書き表される。 なお，x/L は角度[radian]を表す。 この式からディメンション（各項の物理的な意味） を除いて最初の式に変換しても， 振動の本質は失なわない。 なお，最初の式中の a はパラメタであり， 実世界の空気抵抗や重力加速度や振り子の質量と関連する値である。この場合，安定点は， x = n pi (n は整数)になる。 ただし，n が奇数の場合は振り子が上で止まった場合に対応し， 鞍点のため不安定点となり， 本当に安定するのは n が偶数の場合である。 さて，これを解くにあたって次の置き換えを考える： 「2次の微分方程式は，2つの1次の微分方程式を連立させたものである」 ここで速度 v を導入することにより，上の式は次のように書き換えることができる。 v' = x'' = - a * v - sin(x) x' = v さらに，x→x2，v→x1 に置き換えると， x1'=- a * x1 - sin(x2) x2'=x1 となる。この式は，左辺は １ 次の微分だけであり， 右辺は通常の項だけから成る。 これは以下のようにプログラムにできる。 なお，今回は関数内で t を使っていないが， 一般に微分方程式は時間の関数でもあることから， 引数の中には時間 t も入れておく。 ＜myfunc4.m＞ (Download) function xdot=myfunc4(t,x) xdot=zeros(2,1); % 領域確保 a=0.5; xdot(1)=-a*x(1)-sin(x(2)); xdot(2)=x(1); % end of function （Scilab）＜myfunc4.sce＞ (Download) function xdot=myfunc4(t,x) xdot=zeros(2,1); a=0.5; xdot(1)=-a*x(1)-sin(x(2)); xdot(2)=x(1); endfunction さて，まずはこの微分方程式を解くにあたって， まず 変位-速度 平面上での解軌跡を考える。 以下のプログラムを打ち込んで実行してみよう。 この図表から解のふるまいがかなり予想できる。 % prog1e1 - 変位-速度平面上での解軌跡 % 横軸を変位 x ... x2， % 縦軸を速度 v ... x1 とする % x - v 平面を設定します。 x=-9:0.5:9; v=-5:0.5:5; dx=zeros(length(v),length(x),2); for nx=1:length(x) for nv=1:length(v) dx(nv,nx,:)=myfunc4(0,[v(nv) x(nx)]); end end % 各座標での解の動き方を表示します quiver(x, v, dx(:,:,2), dx(:,:,1)) %[xmesh,vmesh]=meshgrid(x, v); %quiver(xmesh, vmesh, dx(:,:,2), dx(:,:,1)) % end of program // prog1e1s (Scilab) // x...x2, v...x1 clear myfunc4 getf('myfunc4.sce'); x=-9:0.5:9; v=-5:0.5:5; dx=zeros(length(x),length(v),2); for nx=1:length(x) for nv=1:length(v) dx(nx ,nv ,:)=myfunc4(0,[v(nv) x(nx)]); end end xbasc(); champ(x, v, dx(:,:,2), dx(:,:,1)) % prog1e1-ver2 - 変数名を変更しただけ % 横軸は変位 x ... x2， % 縦軸は速度 v ... x1. x1=-9:0.5:9; x2=-5:0.5:5; dx=zeros(length(x2),length(x1),2); for nx=1:length(x1) for nv=1:length(x2) dx(nv,nx,:)=myfunc4(0,[x2(nv) x1(nx)]); end end quiver(x1, x2, dx(:,:,2), dx(:,:,1)) %[xmesh,vmesh]=meshgrid(x1, x2); %quiver(xmesh, vmesh, dx(:,:,2), dx(:,:,1)) % end of program // prog1e1s-ver2 (Scilab) // x...x2, v...x1 clear myfunc4 getf('myfunc4.sce'); x1=-9:0.5:9; x2=-5:0.5:5; dx=zeros(length(x1),length(x2),2); for nx=1:length(x1) for nv=1:length(x2) dx(nx ,nv ,:)=myfunc4(0,[x2(nv) x1(nx)]); end end xbasc(); champ(x1, x2, dx(:,:,2), dx(:,:,1)) % prog1e1-ver3 - MATLAB 専用のやり方 % 横軸は変位 x ... x2， % 縦軸は速度 v ... x1. [xmesh,vmesh]=meshgrid(-9:0.5:9, -5:0.5:5); ss=size(xmesh); dx=zeros(ss(1),ss(2),2); for n2=1:ss(2) for n1=1:ss(1) dx(n1,n2,:)=myfunc4(0,[vmesh(n1,n2) xmesh(n1,n2)]); end end quiver(xmesh, vmesh, dx(:,:,2), dx(:,:,1)) % end of program MATLAB と Scilab を比べると， 同じ機能の関数同士でもパラメタの与え方が異なることがある。 ここに示した例では，出力した図は同じであるにもかかわらず， dx というマトリクスのサイズは MATLAB では 21 * 37 * 2 であるのに対し， Scilab では 37 * 21 * 2 である。 それでは，この微分方程式を， 組み込み関数を使わずに解いてみる。 今回は横軸を時刻としてグラフを書く。 % prog1e2 - 振り子の振動 % 計算結果を入れる領域の確保 xx=zeros(501,2); tt=zeros(501,1); % initial condition xx(1,1)=4; % 初期速度 xx(1,2)=0; % 初期変位 tt(1)=0; dt=0.04; % 時間領域: 0.04*500=20(秒) % calculation t=0; for k=1:500 xdot=myfunc4(t,[xx(k,1) xx(k,2)]); xx(k+1,:)= xx(k,:)+dt*xdot'; tt(k+1)=tt(k)+dt; end % display plot(tt,xx(:,1),'b',tt,xx(:,2),'r') % end // prog1e2s (Scilab) // x...x2, v...x1 clear myfunc4 getf('myfunc4.sce'); xx=zeros(501,2); tt=zeros(501,1); xx(1,1)=4; xx(1,2)=0; tt(1)=0; dt=0.04; t=0; for k=1:500 xdot=myfunc4(t,[xx(k,1) xx(k,2)]); xx(k+1,:)= xx(k,:)+dt*xdot'; tt(k+1)=tt(k)+dt; end xbasc(); plot2d(tt,xx(:,1),1); plot2d(tt,xx(:,2),2); MATLAB の組み込み関数 ode23 を使うと， このような微分方程式を簡単に解くことができる。 早速"help ode23"と打ってみよう。 他にも，ode45のように， 目的とインターフェースは同じだがアルゴリズムが微妙に異なる関数が 用意されている。 なお，ここでは，グラフを2回書く。 一枚目は prog1e2 と同じである。 二枚目は横軸を変位，縦軸を速度としたもので， prog1e1 に対応するものである。 一枚目のグラフを書いたら pause で停止するので， (c/r) を （，Scilab では return(c/r) を） 打って，二枚目のグラフに進む。 % prog1e3 - 振り子の振動 % 初期条件は prog1e2 と同じ。 % 時間領域: 0 --- 20 (秒) % 初期速度 4 % 初期変位 0 [tt,xx]=ode23('myfunc4',[0 20], [4 0]); % display plot(tt,xx(:,2),'b',tt,xx(:,1),'r') pause plot(xx(:,2),xx(:,1)) % end // prog1e3s (Scilab) // x...x2, v...x1 clear myfunc4 getf('myfunc4.sce'); tt=0:0.1:20; xx=ode([4 ; 0], 0, tt, myfunc4); xbasc(); plot2d(tt,xx(1,:),1); plot2d(tt,xx(2,:),2); pause; xbasc(); plot2d(xx(2,:),xx(1,:)); （Scilab は MATLAB と違い，ode という関数のみ存在する。 パラメタの与え方も異なる） せっかく解いたので， 最初に示した平面上に解軌跡を表示させてみる。 ここで，hold on というコマンドを使い， 別々に求めた解を，同じグラフ上に表示する。 % prog1e4 - 変位-速度平面上での解軌跡 % 横軸は変位 x ... x2， % 縦軸は速度 v ... x1. [x,v]=meshgrid(-9:0.5:9, -5:0.5:5); ss=size(x); dx=zeros(ss(1),ss(2),2); for n2=1:ss(2) for n1=1:ss(1) dx(n1,n2,:)=myfunc4(0,[v(n1,n2) x(n1,n2)]); end end quiver(x, v, dx(:,:,2), dx(:,:,1)) % 時間領域: 0 --- 20 (秒) % 初期速度 4 % 初期変位 0 [tt,xx]=ode23('myfunc4',[0 20], [4 0]); % グラフの重ね合わせ hold on plot(xx(:,2), xx(:,1),'b') hold off % end // prog1e4s (Scilab) clear myfunc4 getf('myfunc4.sce'); xbasc(); // from prog1e1s x=-9:0.5:9; v=-5:0.5:5; dx=zeros(length(x),length(v),2); for nx=1:length(x) for nv=1:length(v) dx(nx ,nv ,:)=myfunc4(0,[v(nv) x(nx)]); end end champ(x, v, dx(:,:,2), dx(:,:,1)) // from prog1e3s tt=0:0.1:20; xx=ode([4 ; 0], 0, tt, myfunc4); plot2d(xx(2,:),xx(1,:),2); 以上のように，この節では微分方程式の解き方について， 2通りの解き方を示した。 自分で方程式を解く方法は，微分方程式の考え方の習得に有効であった。 一方，ODE23 など MATLAB の関数を使う方法は，極めて少ないプログラミングでも， 実用的な問題解決に有効であることを示した。 ここで示さなかった微分方程式に対しても同様に取り組むことが出来るので参考にして欲しい。

  
  
  

  応用---微分方程式-2 さまざまな振動 （Scilabに違いがあれば）

  先の題材（振り子の運動）を元に， 外力が加わった場合と，パラメトリック振動を取り扱う。 物理的な意味を考えるため，方程式は敢えて物理的なものに遡る。 なお，方程式は次のものに書き換る。 m L x'' = - a L x' - g sin(x) + F(t) 変数 x [rad] は振り子の角度， パラメタ L [m] は振り子の棒の長さ， 記号 ' は時間による微分， パラメタ m [kg] は質量， パラメタ a [kg/s] は空気抵抗， パラメタ g [kg m/s2] は重力加速度， 関数 F(t) [kg m/s2] は外力を表す。 この場合，安定点は， x = n pi (n は整数)になる。 ただし，n が奇数の場合は振り子が上で止まった場合に対応し， 鞍点のため不安定点となり， 本当に安定するのは n が偶数の場合である。先の例題のように， 「2次の微分方程式は，2つの1次の微分方程式を連立させたもの」 と考え，角速度 w を導入するとともに，さらに，x→x2，w→x1 に置き換えると， x'' = w' = x1' = - x1 * a / m - sin(x2) * g / (m * L) + F(t) / (m * L) w = x2' = x1 となる。 外力が加えられた場合 まず関数を用意する。myfunc4 を僅かに手直しするだけである。 関数に引き渡されるパラメタの一つに時刻 t を用意してあるのは， 前出の myfunc4 の時と同様であるが，今回は t を本当に利用している。 それでは早速解いてみる。 次のプログラムは解を求めるものである。 殆んどは前出の myfunc4, prog1e4 と同じであり， 変更点は色違い(赤)で示した。 ただし，コメント行については変更があっても特に示してない。 → ＜myfunc5.m＞ (Download) function xdot=myfunc5(t,x) xdot=zeros(2,1); % 領域確保 F=0.01*sign(sin(0.998*t)); g=1; a=0.1; m=1; L=1; xdot(1)=(-x(1)*a + (-sin(x(2))*g + F)/L )/m; xdot(2)=x(1); % end of function % prog1e5 - 角度-角速度上の解軌跡と，時間-角度，角速度のグラフ % 横軸 は 角度 x ... x2， % 縦軸は角速度 w ... x1. [x,w]=meshgrid(-.15:0.05:.15, -.15:0.05:.15); ss=size(x); dxp=zeros(ss(1),ss(2),2); dxm=zeros(ss(1),ss(2),2); for n2=1:ss(2) for n1=1:ss(1) dxp(n1,n2,:)=myfunc5(0.001,[w(n1,n2) x(n1,n2)]); dxm(n1,n2,:)=myfunc5(-0.001,[w(n1,n2) x(n1,n2)]); end end figure(1); quiver(x, w, dxp(:,:,2), dxp(:,:,1),'r'); hold on; quiver(x, w, dxm(:,:,2), dxm(:,:,1),'g'); hold off % 時間領域: 0 --- 200 (秒) % 初期速度 0, 初期変位 0.01 [tt,xx]=ode23('myfunc5',[0 200], [0 0.01]); % グラフの重ね合わせ hold on plot(xx(:,2), xx(:,1),'b') hold off % 解の時間変化をグラフ化。 figure(2); % 2 枚目のグラフ用紙を用意する plot(tt,xx); % 振動が安定した後の振幅を求める。 % ただし，表示の半分の時刻で振動が安定したと考えた。 len=length(tt); len2=round(len/2); disp('Amplitude = '); disp(max(xx(len2:len,2))-min(xx(len2:len,2))); % end 関数内の外力 F は大きさは時間とともに変わる。 従って，変位-速度平面上のベクトルの方向は時間とともに変わるので注意が必要である。 今回の力は sign (符号) 関数で与えられるため，方形波である。 従って，２通りの力の掛かり方しかないため，ベクトルの図は必要十分の２枚書くことにした。 ２つのベクトルの方向が同一の場合は， 後から書いた矢印によって，最初に書いた矢印の一部が消されているので， しまうので注意しよう。 ここで，上記とは異なる非常に長い周期の外力を加えた場合の解を示す。 二通りの外力に応じた二通り場所に収束しているのがわかる。 myfunc5, prog1e5 からの変更点は次の1箇所である。 myfunc5 3行目 前出の例 F=0.01*sign(sin(0.998*t)); ここでは F=0.01*sign(sin(0.03* (t-1) )); こうした外力を受け，上記のプログラムをそのまま走らせた解を示す。 先ほどよりも大きな軌道を描いている。 これは，振り子独自の振動数にちょうど対応する周期の外力を加えたことで， 大きな振幅が得られたのである。 もしも myfunc5 内の 0.967 という数字を，大きくしても小さくしても， 振幅はこの図よりも小さくなる。 （実際に，0.99 とか 0.9 といった値を当てはめて実行させてみよう。） このとき，角度や角速度の時間変化は，次の図の通りである。 なお，初期値が 0.01 のグラフは角度を，0 のグラフは 角速度を表す。 (最初の図は，収束する様子を示す。後ろの図は，発振が始まる様子を示す。） 基本的にグラフは滑らかな曲線であるが，所々に折れ曲がっているような場所が見られる。 これは，外力が変化したところである。 まとめ：外力があったとき， ・定常状態では外力と振子の周期は等しい ・ちょうどよい周期のとき，振子の振幅は最大になる ・振子の性質（質量，振子の長さ，空気抵抗）に応じて，周期が変わる パラメトリック振動 注意：パラメトリック振動については現在勉強中です。以下の記述は，正しいとは限らない。 パラメトリック振動は，パラメタが時間とともに変化したとき生じる振動である。 振り子の運動に照らし合わせると，時間とともに，例えば振り子の長さ L が変わることに対応する。 振り子とともに長さが変わることは，ブランコを漕ぐ運動の非常に単純なモデル化である。 それでは早速解いてみる。 殆んどは前出の myfunc5, prog1e5 と同じであり， 変更点は色違い(赤)で示した。 ただし，コメント行については変更があっても特に示してない。 ＜myfunc6.m＞ (Download) function xdot=myfunc6(t,x) xdot=zeros(2,1); % 領域確保 F=0; g=1; a=0.1; m=1; L=1-.2*sign(sin(2.04*t)); xdot(1)=(-x(1)*a + (-sin(x(2))*g + F)/L )/m; xdot(2)=x(1); % end of function % prog1e6 - 角度-角速度上の解軌跡と，時間-角度，角速度のグラフ % 横軸 は 角度 x ... x2， % 縦軸は角速度 w ... x1. [x,w]=meshgrid(-.8:0.2:.8, -.8:0.2:.8); ss=size(x); dxp=zeros(ss(1),ss(2),2); dxm=zeros(ss(1),ss(2),2); for n2=1:ss(2) for n1=1:ss(1) dxp(n1,n2,:)=myfunc6(0.001,[w(n1,n2) x(n1,n2)]); dxm(n1,n2,:)=myfunc6(-0.001,[w(n1,n2) x(n1,n2)]); end end figure(1); quiver(x, w, dxp(:,:,2), dxp(:,:,1),'r'); hold on; quiver(x, w, dxm(:,:,2), dxm(:,:,1),'g'); hold off % 時間領域: 0 --- 500 (秒) % 初期速度 0, 初期変位 0.01 [tt,xx]=ode23('myfunc6',[0 500], [0 0.01]); % グラフの重ね合わせ hold on plot(xx(:,2), xx(:,1),'b') hold off figure(2); plot(tt,xx); % 振動が安定した後の振幅を求める。 % ただし，表示の半分の時刻で振動が安定したと考えた。 len=length(tt); len2=round(len/2); disp('Amplitude = '); disp(max(xx(len2:len,2))-min(xx(len2:len,2))); % end 時間とともにゆれが大きくなり，最後には振幅が安定する。 まとめ：ブランコの長さが変わるというパラメトリック振動では， ・定常状態ではブランコの伸び縮みの周期は，ブランコの揺れ周期の半分である ・このブランコの場合，伸び縮みの関数として， L=1-.2*sign(sin(2.04*t)); L=1-.2*sign(sin(1.96*t)); くらいの範囲がちょうどよい周期である。 ブランコは振れが大きいほど周期が長くなる性質があり， 励起の周期が長い（カッコ内の数字が 1.96 に近い）程，最終的な揺れの振幅は大きいが， 安定点に到達するまでの時間が長い。 一方，励起の周期が短い（カッコ内の数字が 2.04 に近い）程，最終的な振幅は小さいが， 安定点に到達するまでの時間は短い。 ここで示した数字よりも長すぎたり短すぎる励起の周期の場合は， ブランコはあまりよく振れてくれないようである。 ただし，初期条件によって，これらの性質は変わる可能性がある。 ・ブランコの性質（質量，振子の長さ，空気抵抗）に応じて，周期が変わる 以上のように，この節では振動が継続する場合について， 2通りの例題を示した。 ここでは表示をきれいにするためのプログラムを入れたが， そういうことが不要であれば，極めて少ないプログラミングで， 実用的な問題解決に有効である。

  
  

  ---

  応用1f --- s を使って表記された特性のグラフ化

  回路解析をしていると，伝達関数を s を使って表すことがよく行われます。 s は微分差要素です。ここではこれを使った関数をグラフ化してみましょう。
  

  応用---s を使って表記された特性のグラフ化 （Scilabに違いがあれば）

  RC 回路の時間応答 (Cの電圧) 次の回路は，応用1aの節でも取り扱ったものである。 入力 Vin を電圧 V1，出力 Vc をコンデンサにかかる電圧とすると， 出力と入力の関係式は， Vc = 1 1 + τ s Vin で表わすことができる。ただし，τ=CR である。 さて，s は微分の演算子である。 このことから，微分方程式を解く事が，そのまま回路特性を求めることになる。 これ以降特性を求めるが，特に断りがない限り， 「ステップ関数の入力電圧に対する出力」を求めることにする。 ということは，時刻 t=0 以降は Vin=1 (固定) ということである。 さて，これを解くにあたって，最初に次の置き換えを行う： 出力：Vc → y 続いて，方程式を次のように書き直す。 ただし，記号のダッシュ( ' )は時間による微分を表すものとする。 経過 ・y + τ y' = 1 結論 ・y' = (1 - y)/τ ここまで置き換えれば，前節の手順をそのまま使って解く事ができる。 ＜myfunc7.m＞ (Download) function ydot=myfunc7(t,y) ydot=zeros(1,1); % 領域確保 tau=1000 * 1.e-6; ydot(1)=(1-y(1))/tau; % end of function （Scilab）＜myfunc7.sce＞ (Download) function ydot=myfunc7(t,y) ydot=zeros(1,1); tau=1000 * 1.e-6; ydot(1)=(1-y(1))/tau; endfunction それでは，早速解いて見よう。 % prog1f1 - RC回路の応答 % 時間領域: 0 --- 10 (ms) % 初期電圧 0 [tt,yy]=ode23('myfunc7',[0 10.e-3], [0]); vc=yy(:,1); plot(tt,vc) // prog1f1s (Scilab) clear myfunc7 getf('myfunc7.sce'); tt=0:1.e-4:10.e-3; yy=ode([0], 0, tt, myfunc7); vc=yy(1,:); xbasc(); plot2d(tt,vc,1); RC 回路の時間応答 (Rの電圧) 同じ回路を用いて， 抵抗にかかる電圧 Vr を出力とした場合についてグラフ化しよう。 式は， Vr = τ ｓ 1 + τ ｓ Vin = { 1 - 1 1 + τ ｓ } Vin で表わすことができる。ただし，τ=1/(CR) である。 二つ目の等号から，Vr ← Vin - Vc であることがわかる。 従って，prog1f1 の続きとして， 次の命令を続ければ，Vr をグラフ化することができる。 % prog1f1 の続き vr=1-vc; plot(tt,vc,tt,vr) // prog1f1s の続き(Scilab) vr=1-vc; plot2d(tt,vr,2); 2次 の LPF の 時間応答 （LPF とは Low Pass Filter のこと） 続いて，応用1a の節でも取り扱った 2次の回路について特性をグラフ化する。 入力 Vin を電圧 V1，出力 VC2 をコンデンサ C2 にかかる電圧とすると，VC2 は， VC2 = 1 1 + a * s + b * s^2 Vin で表わすことができる。ただし， a = C2R2 + C2R1 + C1R1 b = C1C2 R1R2 である。また，ステップ入力を考えるため以降 Vin=0 とする。 この特性は，低い周波数では伝達率が 1 (=100%) だが， 高い周波数では伝達率が小さくなることから，LPF の特性と言える。 先の時と同様に， 「s は微分の演算子」ということを考慮し， 最初に次の置き換えを行う：VC2 → y 続いて，方程式を次のように書き直す。 ただし，記号のダッシュ( ' )は時間による微分を表すものとする。 経過 ・y + a s y + b s^2 y = 1 経過2 ・y2 = s y ... ＜y2の導入＞ ・y + a y2 + b s y2 = 1 経過3 ・y2 = y' ・y + a y2 + b y2' = 1 結論 ・y' = y2 ・y2' = 1/b - y/b - y2*a/b ここまで置き換えれば，前節や直前の手順をそのまま使って解く事ができる。 ＜myfunc8.m＞ (Download) function ydot=myfunc8(t,y) ydot=zeros(2,1); % 領域確保 a=1.e-6*100.e3 + 1.e-6*1000 + 0.1e-6*1000; b=1.e-6*0.1e-6 * 100e3*1000; ydot(1)=y(2); ydot(2)=(1 - y(1) - y(2)*a)/b; % end of function （Scilab）開発中 % prog1f2 - RC回路の応答 % 時間領域: 0 --- 10 (ms) % 初期電圧 0 [tt,yy]=ode23('myfunc8',[0 5.e-3], [0 0]); vc2=yy(:,1); plot(tt,vc2) (Scilab 開発中) ちなみに，VC2の過渡応答を解釈するには， ２つの時定数を考えればよい。 一つは，R1，C2 から成る回路による 時定数であり，τはおよそ 0.1 ms である。 もう一つは，R2，C2 から成る回路による 時定数であり，τはおよそ 100 ms である。 従って，グラフを書く際の範囲の選び方として， 最初の 0.5ms 程度とすれば， R1，C1 による効果がわかるし， 500ms 程度とすれば， R2，C2 による効果がわかる。 2次の BPF の時間応答 （BPF とは Band Pass Filter のこと） 続いて，出力として，同じ回路の抵抗 R2 として計算する。 VR2 は， VR2 = s C2 R2 1 + a * s + b * s^2 ・Vin = s C2 R2 ・VC2 で表わすことができる。 この特性は，低い周波数と高い周波数で伝達率が小さくなるが， 最大の伝達率はその中間の適当な周波数であり，正に BPF の特性と言える。 先のprog1f1 にて VC2 は既にグラフ化した。 また，先立つ prog1f1 の「経過2」の計算から， s * VC2 は y2 即ち，yy(:,2) として求められる。 従って，prog1f2 の続きとして， 次の命令を続ければ，VR2 をグラフ化することができる。 % prog1f2 の続き vr2=100.e3 * 1.e-6 * yy(:,2); plot(tt,vc2,tt,vr2) // prog1f2s の続き(Scilab)＜開発中＞ ちなみに，VR2 の過渡応答を解釈するのに，R1とC2 から得られるおよそ0.1msの時定数と，R2とC2 から得られるおよそ 100 msの時定数を意識する必要があるのは，VC2 の過渡応答を見るときと同様である。 2次の HPF の時間応答 （HPF とは High Pass Filter のこと） 続いて，出力として，次のような特性を考えてみる。 Vout ∝ s^2 ・VC2 = s^2 1 + a * s + b * s^2 ・Vin この特性は，分子に s^2 に比例する項がある。 この項は，計算によって求めた yy( ) を使って求めることはできない。 HPF の特性を得るためには，1 から，s の 0 乗や s の 1 乗の項を引くことである。 従って，prog1f2 の続きとして， 次の命令を続ければ，Vhpf をグラフ化することができる。 % prog1f2 の続き a = 1.e-6*100.e3 + 1.e-6*1000 + 0.1e-6*1000; vhpf = 1 - yy(:,1) - a * yy(:,2); plot(tt,vc2,'K',tt,vr2,'B',tt,vhpf,'R') % blacK, Blue, Red // prog1f2s の続き(Scilab)＜開発中＞ ※ ode23 による計算のためのパラメタとして，時間領域として，0.5 ms と 500 ms それぞれを設定してみよう。 2次 系の時間応答 おさらい おさらいとして，2次の系の時間応答をまとめる。 ここでは，L-R-C　直列回路を仮定する。 先ずは LPF の特性として，コンデンサの電圧を考える。 VC = 1 1 + s * a1 + s^2 * a2 ・Vin ただし， a1 = C * R , a2 = C * L である。 「s は微分の演算子」ということを考慮して書き直す。 ただし，VC → y とし，記号( ' )は時間微分とする。 経過 ・y2 = s y ... ＜y2の導入＞ ・y + a1 y2 + a2 y2' = 1 結論 ・y' = y2 ・y2' = 1/a2 - y/a2 - y2*a1/a2 定数は，次の通りとする。 C=1.e-9[F] R=1000[Ω] L=0.001[H] ＜myfunc9.m＞ (Download) function ydot=myfunc9(t,y) ydot=zeros(2,1); % 領域確保 a1=1.e-9 * 1000; a2=1.e-9 * 0.001; ydot(1)=y(2); ydot(2)=(1 - y(1) - y(2)*a1)/a2; % end of function （Scilab）開発中 % prog1f3 - RLC回路の応答 % 時間領域: 0 --- 30 (us) % 初期電圧 0 [tt,yy]=ode23('myfunc9',[0 30.e-6], [0 0]); vc=yy(:,1); plot(tt,vc,'K'); % K---blacK (Scilab 開発中) 続いて，出力を同じ回路の抵抗 R として計算する。 VR は， VR = s C R ・ VC で表わすことができる。 先立つ prog1f3 の「経過」の計算から， s * VC は y2 即ち，yy(:,2) として求められる。 % prog1f3 の続き hold on vr=1.e-9 * 1000 * yy(:,2); plot(tt,vr,'B'); % B---Blue // prog1f2s の続き(Scilab)＜開発中＞ 続いて，出力を同じ回路のコイル L として計算する。 この項を求めるには，1 (=Vin) から，抵抗やコンデンサの電圧を引くことである。 % prog1f3 の続き vL = 1 - vc - vr; plot(tt,vL,'R'); % R---Red hold off // prog1f3s の続き(Scilab)＜開発中＞ ※ 確認しよう ： R の値を変えると，減衰率が変化する。 CまたはLを変えると振動の周期が変化する。 ※ Vin がランプ波（１次式で時間と共に増加する波形）だった場合， は，本来 Vin が収まる場所として 定数×t を挿入すれば良い。 他にも正弦波など様々な入力信号が考えられるが， その部分を適切な関数にすれば良い。 ＜myfunc7 〜 myfunc9 の改良＞ （ランプ波 Vin=30000*t として） function ydot=myfunc7(t,y) ydot=zeros(1,1); % 領域確保 tau=1/(1000 * 1.e-6); ydot(1)=(30000*t - y(1))*tau; % end of function function ydot=myfunc8(t,y) ydot=zeros(2,1); % 領域確保 a=1.e-6*100.e3 + 1.e-6*1000 + 0.1e-6*1000; b=1.e-6*0.1e-6 * 100e3*1000; ydot(1)=y(2); ydot(2)=(30000*t - y(1) - y(2)*a)/b; % end of function function ydot=myfunc9(t,y) ydot=zeros(2,1); % 領域確保 a1=1.e-9 * 1000; a2=1.e-9 * 0.001; ydot(1)=y(2); ydot(2)=(30000*t - y(1) - y(2)*a1)/a2; % end of function （Scilab）開発中

  
  

  ---

  応用1g --- 偏微分方程式（半導体中のキャリアの運動）

  偏微分方程式を解くことだってできます。

  応用---偏微分方程式 （Scilabに違いがあれば）

  ここでは，半導体中のキャリアの運動について，MATLAB を利用して解きます。 ただし，数式の理解のために MATLAB を使うのであって， 各項の大きさについては特に物理的な根拠はありません。 方程式：n型半導体内の少数キャリア（正孔）の運動を表す式 （ここで， ' とは，x による微分を表す） dp/dt = D * p'' - u * E * p' - u * p * E' + G - (p - p0)/T ただし， p は正孔濃度 p0 は熱平衡時の正孔濃度 D は拡散定数 u は移動度 E は電界 G は生成項 T は再結合時間 まずは，このプログラムが動くように設定し， 一度走らせて見よう。 解説については，その後に行う。 % prog1g1 - 半導体中の正孔の振舞 xx=100; x=1:13*xx; lenx=length(x); p0=1; p=ones(1,lenx)*p0; % 100 (1 line) p(6*xx:8*xx)=p0*(1 + ... min([ones(1,2*xx+1) ; ... linspace(0,4,2*xx+1) ; ... linspace(4,0,2*xx+1)])/2); % 110 (2 lines) p(1*xx:2*xx)=p0+2*p0; p(2*xx:3*xx)=0.5*p0; ee=zeros(1,lenx); % 120 (3 lines) ee(5*xx:9*xx)=1; ee(4*xx:5*xx)=linspace(0, 1, xx+1); ee(9*xx:10*xx)=linspace(1, 0, xx+1); g=zeros(1,lenx); % 130 (1 line) g(11*xx:12*xx)=1; d=0.2; hold off plot(x,p) axis([0,lenx,0,7.5]); hold on; plot(x,4.5+ee,'g'); ee=0.02*ee; plot(x,6+g,'y'); g=0.0005*g; % 140 (2 lines) for aa=1:3 for a=1:100 p(1)=0; p(lenx)=0; pp=p*(1 - 2*d); pp=pp + [0 p(1:lenx-1)].*(d + [0 ee(1:lenx-1)]); pp=pp + [p(2:lenx) 0].*(d - [ee(2:lenx) 0]); % one of these two % 210 (2 lines) pp(2)=pp(2)+pp(1); % pp(2)=0; % one of these two % 220 (2 lines) pp(lenx-1)=pp(lenx-1)+pp(lenx); % pp(lenx-1)=0; p=pp; % 230 (2 lines) p=p+g; p=p-0.00015*(p-p0); end pause(0.5); plot(x,p); end plot(x,p,'r'); disp('--end--'); % end of program // prog1g1s (Scilab) xx=100; x=1:13*xx; lenx=length(x); p0=1; p=ones(1,lenx)*p0; // 100 (1 line) for a=6*xx:8*xx p(a)=p0*(1+min([0.5 (a-6*xx)/xx (8*xx-a)/xx])); end // 110 (2 lines) p(1*xx:2*xx)=p0+2*p0; p(2*xx:3*xx)=0.5*p0; ee=zeros(1,lenx); // 120 (3 lines) ee(5*xx:9*xx)=1; ee(4*xx:5*xx)=linspace(0, 1, xx+1); ee(9*xx:10*xx)=linspace(1, 0, xx+1); g=zeros(1,lenx); // 130 (1 line) g(11*xx:12*xx)=1; d=0.2; xbasc(); plot2d(x,p); plot2d(x,4.5+ee); ee=0.02*ee; plot2d(x,6+g); g=0.0005*g; // 140 (2 lines) for aa=1:3 for a=1:100 p(1)=0; p(lenx)=0; pp=p*(1 - 2*d); pp=pp+[0 p(1:lenx-1)].*(d+[0 ee(1:lenx-1)]); pp=pp+[p(2:lenx) 0].*(d-[ee(2:lenx) 0]); // 210 (2 lines) pp(2)=pp(2)+pp(1); // pp(2)=0; // 220 (2 lines) pp(lenx-1)=pp(lenx-1)+pp(lenx); // pp(lenx-1)=0; p=pp; p=p+g; p=p-0.00015*(p-p0); end plot2d(x,p); end plot2d(x,p,2); disp('--end--'); // end of program MATLAB → Scilab の際に， % コメント → // コメント， plot → plot2d， hold 関数 → xbasc という変更がなされた。しかし，殆んどは同一のプログラムである。 それでは，この微分方程式を微妙に変えてみよう。 注意：この計算の中で，ドリフト電流を計算する際には， 拡散項は省くことができない。 % 100 の次から 3 行は，電界中でのキャリア濃度の変化を表す。 この 3 行を消すと，その部分のキャリアを消せる。 % 110 の次から 2 行は，電界部の左側のキャリア濃度の変化を表す。 この 3 行を消すと，その部分のキャリアを消せる。 % 120 の次から 3 行は，電界部を表す。 この 3 行を消すと，電界をなくすことができる。 % 130 の次の 1 行は，生成項を持つ領域を表す。 この行を消すと，生成項を消せる。 % 140 の次の 2 行は，繰返しの回数を指定する。 aa に関する繰返しは，プロットする回数を指定する。 この繰返し回数を増やすと，長時間にわたるシミュレーションができる。 a に関する繰返し回数は，プロットする間隔を指定する。 % 210 の次から 2 行は， 系の左端での拡散現象に関する境界条件を表す。 2 行のどれもがコメントアウトしてもしなくても良い。 どちらもコメントアウトしたときは，ちょうど窓が開いているように， キャリアが自然に外に漏れてゆくことを示す。 1 行目のみ生かして2行目をコメントアウトしたときは， ちょうど窓が閉じているように，キャリアが外に漏れないことを示す。 2 行目の物理的な意味は，強制的にその場の濃度が決まることを意味する。 プログラム例では 0 に固定しているが，ほかの値を設定しても良い。 2 行とも生かしたときは，2行目だけを生かした場合と同様。 % 220 の次から 2 行は，% 210 と同様であるが， 系の右端での拡散現象に関する境界条件を表す。 % 230 の次から 2 行は， それぞれ生成項と，生滅項の計算部である。 これらの行をコメントアウトすれば， その物理現象のシミュレーションを省くことができる。