back サーバ選択 (基本ミラー)
Last Update : '16.4.24 「応用1a」内のミスプリ1個所の修正と,「外部リンク 応用1a-2016」を追加; '07.5.8言い回しの修正; '07.2.15 例題の番号を振りなおした; '07.1.31 Filter を Scilab 化できた; '07.1.24 = ミスプリ訂正(Change Directry のメニューが「edit」だったのを「file」に正した 等); '07.1.23 = 演算C をスキップしながら,応用1b まで焼きなおした; '07.1.18 = 「matlabつかいませんか」をもとに,「Scilabつかいませんか」をつくりはじめる; '05/12/20 = 「s で表現された特性の応答」に関しパラメタ τ の使い方を修正 ; --- 略 --- ; '98/12/11 = 微分方程式の項を書き加えた。
Scilab つかいませんか

--- Scilab を使うための ABC ---
★注意:'07.5.8現在,
「matlabつかいませんか」から
「Scilab つかいませんか」への
移行作業は残り約 25% です★
Link
目次
  • Scilab Overview
  • 実行例
  • Scilab 環境の説明 (ディレクトリ,ファイル,sceファイル)
  • 演算
    • 演算a - 電卓として使う,関数一覧
    • 演算b - 1次元配列,配列のグラフ化
    • 演算c - 文字列
    • 演算d - 2次元配列,配列のグラフ化
    • 演算e - マトリクス
  • プログラミング(sce-ファイル)
  • 応用 1
    • 応用1a - 交流回路の電圧波形(計算とグラフ化)
         【外部リンク】応用1a-2016 PBL2016の開発に使う例
    • 応用1b - 実験データの整理(データと,回帰線)
    • 応用1c - Wave ファイル と 信号処理 (FFT,フィルタ) 
    • 応用1d - 数値積分
    • 応用1e - 微分方程式 (振子の運動, 外力があったとき, パラメトリック振動) (←半分matlab残る
    • 応用1f - s を使って表記された特性のグラフ化 (←まだmatlabのまま
    • 応用1g - 偏微分方程式 (半導体中の少数キャリアの運動)
  • 応用 2 (link) ( Sorry! Staff Only)
    • 電子回路や電気回路の特性を Scilab で数値計算する (←まだmatlabのまま
  • 応用 3 (link)
    • 電子回路や電気回路の特性を Symbolic に求める (←まだmatlabのまま
  • 応用 4 ( Sorry! Staff Only )
    • ラプラス変換を使った解析 (←まだmatlabのまま


Scilab Overview

私が Scilab について話すには,まず MATLAB に触れなくてはなりません。

MATLABは,FORTRANを基に, マトリクス演算を強力に扱えるようにしたコンピュータ言語です。
MATLABは,次のような特徴を持っています:
[インタープリタ型である]そのため, ちょうど電卓を使うように, コマンドラインから1行ずつコマンドを打ち込んで, その都度実行することが可能です。 この機能のお陰で,MATLABを使って問題解決に取り組んだユーザは, プログラムを作らなくても,かなりの問題をとくことができます。
[MATLABはマトリクスや複素数を扱う数値計算が大得意です] もともとのMATLABという名前にもその特徴が現れています。 マトリクスの演算は,普通の足し算や掛け算の記号を使うだけでOKです。
[MATLABはコンピュータ言語としても使用可能です] ループ・判断・サブルーチンを含む一連のプログラムによって 複雑な処理もこなすことができます。 しかもインタープリタ型なので,デバッグは非常に簡単です。


さて,ようやく Scilab の説明を始めましょう。 Scilab は,MATLABとほぼ互換のフリーのソフトのひとつです。 互換の特長から,MATLABの利点はそのままScilabの利点になります。 また,フリーですから,自分のパソコンに自由にインストールできますし, 学校の多くのパソコン上にもインストール済みです。 しいて欠点をあげれば,日本語の解説が不足していることだと思います。

このページは,望月が「自分が使い方を忘れないため」に作るページです。 また,私の研究室の学生に Scilab を学んでもらうために作るページでもあります。 従って,例題などは私の研究に近いところでだけ作っています。 ですから,他の部署の方々には殆ど役立たない可能性があります。

なお,このページを作るにあたって主に参考にした情報は次のとおりです。 また,関連するリンクは,目次のところにまとめました。

起動

この説明は Scilab 4.1 を前提にします。

インストールするには, Scilab 日本語ページ からダウンロードしてください。 普通にインストールすれば,特に問題なくインストールできると思います。
なお,沼津高専の総合情報センターでは,既にインストール済みですから, 直ぐに使い始めることができます。

「スタート」ボタンから Scilab を探して実行させると,次のようになります。
◎Scilab Command Window [_][□][×]
File Edit Preference Control Editor Applications ?
[ ] [ ] [ ] [ ] [ ] [ ] . . .
        ___________________________________________
                         Scilab-4.1

                  Copyright (c) 1989-2006
              Consortium Scilab (INRIA, ENPC)
        ___________________________________________


Startup execution:
  loading initial environment

--> 





File メニューからは,次のコマンドを実行できます:
new Scilab , exec... , open... , load... , save... , change directory , get current directry , print setup , print , exit

Edit メニューからは,次のコマンドを実行できます:
select all , copy , paste , empty clipboard , history

Preferences メニューからは,次のコマンドを実行できます:
languadge , colors , toolbar , files association , choose font , clear histroy , clear command window , console

Control メニューからは,次のコマンドを実行できます:
resume , abort , interrupt

Editor メニューからは,専用のテキストエディタを実行できます:

Applications メニューからは,次のコマンドを実行できます:
scicos , edit graph , m2sci , browser valiables

? (help) メニューからは,次のコマンドを実行できます:
Scilab help , configure , Scilab demos , web links , aboart



実行例

ここに,Scilab 4.1 を実際に使って簡単な計算をした例を示します。
図の後に,各部を説明します。
◎Scilab [_][□][×]
File Edit Preferences Control Editor Applications ?
--> 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)
  • ここでは,実行画面を示します。
  • 本来はもっと空白行が多いですが,ここでは,空白行の一部は詰めました
  • ここでは入力はこの色で示します。(本来はモノクロ画面)
  • (括弧)はコメント用に私が書き加えたものです
  • (1) : 対話式に入力します。
    普通のコンピュータ言語(BASIC, C など)のように式を書けます。 なお,BASICと同様,定義しない変数を使えます。
    直ぐに計算し,結果を表示します。
    ";"や","で区切れば1行に複数の式も書けます。 (";"区切りは,結果表示無し)
  • (2) : この処理系の特徴である,配列の一括処理の例です。
    0 から 3 までの配列を定義し,要素それぞれに 30 を掛けます。 計算結果も,要素それぞれを示します。
  • (3) : %pi (=π=3.14・・・) / 180 の掛け算も,sin という計算も, 配列全体に一括処理できます。
    30度刻みで sin(0度)から,sin(90度)まで計算出来ました。
  • (4) : 入力プロンプトが出て,次の入力を待っています。 以下,同様に処理を続けられます。


環境の説明(ディレクトリ,ファイル,sce-ファイル)

演算結果をファイルに保存することができます。
現在使用中のディレクトリはカレントディレクトリです。
テキスト形式のファイルで,プログラムを作ることができます。
(Scilab では,.プログラムを sce ファイルに,関数を .sci ファイル に格納するようです)
代表的には以下のようなコマンドを覚えておけば大丈夫。
環境
基本的な操作:
ヘルプ プロンプトからキー入力 -->help
「-->」というプロンプトが出ているとき,「help」とエンターキーを打つ。以下同様。
メニューからの操作  「?」>「help」
メニューの「?」を選択し,続いて「help」を選択する。以下同様。
※ Scilab のヘルプメッセージは英語です。
あいまい検索調査中
MATLABでは「lookcor keyword」にてkeyword に関連するコマンドを調べる
マニュアル プロンプトからキー入力 -->help  <keyword>
(keyword に関するマニュアルを表示)
プロンプトからキー入力の例 -->help sin
ヒストリ プロンプトにてカーソルキー を押せば,前に入力した行が復帰。何行でも戻れるし,再編集できる
メニューからの操作 「Edit」>「History」 
デモ メニューからの操作  「?」>「Scilab Demo」
終了 プロンプトからキー入力 -->exit
メニューからの操作  「File>「exit」
カレントディレクトリ(作業領域)について:
ファイルのリスト プロンプトからキー入力 -->dir                 (コンパクトな形式
プロンプトからキー入力 -->unix_x("dir")   (詳細な形式
※ Scilab の unix コマンドは, OS (含 Win.) のコマンド実行。 _w は,出力を Scilab ウインドウ内に設定
パス表示 プロンプトからキー入力 -->pwd
メニューからの操作  「file」>「Get Current Directory」
パス変更 プロンプトからキー入力 -->chdir <path>   (カレントディレクトリを,path に設定する
メニューからの操作 「file」>「Change Directory」
エディタについて:
エディタ起動 メニューからの操作  「Editor」
エディタ起動後に,操作すべきファイルを指定したり,セーブする
なお,Scilab では拡張子は明示的に指定しなくてはならない。(通常のWindowsソフトなら,例えばテキストエディタなら.txtなどのような拡張子を勝手に付けてくれた) 拡張子を特に指定しない場合は,拡張子のないファイルが出来てしまう。ただし,そういったファイル名を指定すればプログラムは動いてしまう。
なお,標準は *.sce という拡張子なので,拡張子にいちいち .sce としてつけること。関数の場合 *.sci という拡張子にするのが基本。


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

最も簡単な応用方法は,電卓として利用することです。
倍精度型が標準ですから,高い精度の計算が保証されます。
変数を使えは,計算結果を一時的に記憶しておくことができます。 その場合,「変数 = 数式」 という書き方をするのは,他の言語と全く同じです。 この「変数 = 数式」という形式を代入文と呼びます。 前もって定義すること無しに変数を使えます。 なお,一時的というのは, 終了したり,clear コマンドを実行したり,同じ変数に対して改めて代入するまでという意味です。
代表的には以下のようなコマンドや,使い方を覚えておけば大丈夫。
演算---電卓として使う,関数一覧
演算と代入文
ans
入力例  -->1+2  (答え(今回は3)を表示するとともに,値はデフォルト変数「ans」に記録される
入力例  -->ans+4  (ansは,次の計算に使える。今回の答えは7
入力例  -->y=2+3  (イコールから右が計算(評価)され,値は自動的に用意された変数「y」に記録(代入)される
入力例  -->a=y+3  (変数は新たな計算に使うことができる
入力例  -->y=y+1  (数学ではありえない数式だがscilab(を始めとする計算機言語)では許される。この場合,先ずイコールから右が計算(評価)され,値は変数「y」に記録(代入)される。計算の結果,yの値は,さっきまでのyの値よりも1増える
変数のみ + (C/R) 入力例  -->y  (yの内容を表示する
";" 入力例  -->z=y+5;  (最後のセミコロンは「表示無し」の意味
";"と"," 入力例  -->a=y+1; b=y+2, c=y+3   (";"と","は,1行に複数のコマンドを書く場合の区切記号として使われる。この場合";"で区切れば表示無しであり","で区切れば表示する
四則演算とべき乗 a + b, a - b, a * b, a / b (or b \ a), a ^ b
inf 1/0 のこと:無限大
NaN または nan Not a Number : 数字じゃない
円周率 %pi と表記
虚数 %i と表記
(%i ^ 2 = -1) 例:"a=1 - 2*%i"
変数の保管 入力例-->save mywork.mat
mywork.mat というファイルに使用中の変数全体を保管する
Scilabは拡張子を補わないので, .mat ファイルに保存するにはこの例の通りに拡張子付きで mywork.mat と打ち込む
変数の復帰 入力例-->load mywork.mat
mywork.mat というファイルから,変数を復帰する
Scilabは拡張子を補わないので, .mat ファイルから復帰するにはこの例の通りに拡張子付きで mywork.mat と打ち込む
変数名確認 入力例-->who簡略版
入力例-->whos詳細版
変数の破棄 入力例-->clearすべての変数を破棄
入力例-->clear aa という変数を破棄
宣言 C や Pascal と異なり,宣言しない変数を使うことができる。 変数名は,1文字目がアルファベットで, 2文字目以降は数字とアンダースコア'_'も許される。
A と a 大文字と小文字を区別する。A と a は別物であり,同時に共存できる。
Scilabの古いバージョンでは大文字と小文字を区別しないものもあった
数値表示の形式 入力例-->format('v',10)formatコマンドは今後の数値の表示形式を標準とする
入力例-->format('e',15)固定小数点15桁表示
入力例-->format('v',15)固定or浮動小数点15桁表
関数 sin, sqrt など。y=sin(x) など,普通の関数は普通にあります。 ("help sin" など打ち込んで確かめよう)
殆んどの関数も複素演算できるので,電気回路理論と整合が良い。 例:y = real(v * exp(j*2*%pi*f*t))


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

良いところは,配列を簡単に扱える点です。
配列を使って代入文を書くと,自動的に"要素ごとの演算"をやってくれます。
以下にその利点を例を示します。
Scilab BASIC
(対応した書き方)
BASIC
(BASICに最適化)
>> x=linspace(0,%pi,101);

>> y=sin(x);

>> plot2d(x,y);
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
上記の例では Scilab ではずいぶんと簡潔に書けています。 しかも標準のサブルーチン plot のおかげで, 簡単にグラフが描けてしまいます。 ここで BASIC の名誉のために言っておきますが, Scilab と BASIC のどちらも, もう一方に対して劣っているわけではありません。 両方の言語ともに適した分野が有り, たまたま上記のような例は, Scilab が最も得意としたということです。 すなわち,配列を使い,その結果を効率良く表示するという分野です。
ところで,Scilab では,1次元配列を扱うとき, 関数の引数等でちょっと不思議なところがあります。(例 : size 関数) これは,Scilab がもともと MATRIX (2次元配列) を扱うことを前提としているからでしょう。
代表的には以下のようなコマンドや,使い方を覚えておけば大丈夫。
演算---1次元配列,配列のグラフ化 (MATLABに違いがあれば)
x(3) x という配列の3番目の要素
x=[6 9 4] x(1)=6, x(2)=9, x(3)=4 という配列を作る
x(2)=7
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次元のグラフ
scf(n) これ以降のグラフィック関係の表示先として,n という番号のグラフ表示ウインドウにします。デフォルトで scf(0) が指定済み。 
もしも,scf(n) として初めてのウインドウを指定した場合は,まっさらなウインドウが自動的に用意されますが, 2回目以降だった場合は,前のデータをそのまま残した状態です。 もしも新しく書き直したいときは clf 命令を使います。 
clf(n) n という番号のついたグラフ表示ウインドウを,まっさらな状態にします。
もしも clr() という風に引数を指定しない場合,直前の scf(n) で指定したウインドウに作用します。
plot(x,y) 業界標準の MATLAB の標準的な2次元プロット関数に倣って作られた関数です。MATLAB 上の plot 関数とほぼ同様です。
なお,plot 関数はリニア専用です。ログの表示のためには本来のMATLABでは,loglogやsemilogという関数を使いますが,それらについては Scilab は用意してありません。
Scilab の標準的な2次元プロットは plot2d です。この関数は,対数表示の機能まで備えた優れものです。
plot2d(x,y) Scilab 本来の2次元グラフ表示関数。
可能:plto2d(y)   (xとして [1:length(y)] を仮定)

可能:plot2d(x,y,logflag="nn") は,plot2d(x,y,logflag="nn")  と同じです。
可能:plot2d(x,y,logflag="nl") は,縦軸(y軸)をlogにします。
可能:plot2d(x,y,logflag="ln") は,横軸(x軸)をlogにします。
可能:plot2d(x,y,logflag="ll") は,全対数グラフにします。

可能:plot2d(x,y,style=1) は,色を黒にします。数値が色を表します。
0点線,1黒,2青,3緑,4水,5赤,6紫,7黄
これ以外の数値も指定できます。
負の値は,プロットです。

可能:plot2d(x,y,logflag="ll",style=2) という風に複数指定できます。
同じ座標に3つのグラフを表示 最初に次の操作により 3 つの関数を定義する。 これらを一つグラフに書いてみよう。
x=0:0.1:7;
y1=sin(1*x);
y2=sin(2*x);
y3=sin(3*x);
Scilab では,clf により画面を消去する。 plot は簡単な機能しか持たないようであるが, plot2d はそこそこの能力を持っているので plot2d を使う。 plot2d は前の画面を消さずにグラフを書く。 (初めてグラフを書くときは不要だが,通常は clf による画面消去が必要) plot2d の style パラメタは色を指定するが, 負の値を設定するとプロットになる。
なお,plot2d の拡張形は plot2d(x,y,logflag='nn',style=n) となる。
 'nn' は両軸とも通常目盛りを意味する。 デフォルト設定でもあり指定無き場合は'nn'とみなされる。 代わりに 'ln' とすれば横軸だけが Log であり, 'nl' なら縦軸だけが Log であり, 'll' なら両対数グラフとなる。 大文字による指定('NL'など)は禁止である
clf()
plot2d(x,y1,style=1)
plot2d(x,y2,style=2)
plot2d(x,y3,style=3)

これ以降,配列専用の関数-1 (1次元でも良く使うもの)
max, min, sort など max(a) とは a という配列の最大値を求めます
(MATLAB と 微妙な違いあり。焼き直しの際は注意が必要)
length(a) a の要素の数をもとめる。もしも m行×1列 または 1行×m列 のベクトルならば答えは m である。もしも m行×n列の行列なら, m×n が返ってきます。
MATLAB と 微妙な違いあり。焼き直しの際は注意が必要
size(a) a という配列の行と列を求めます。Scilab は標準的な数値を2次元配列だと仮定してるので,もしもスカラーでも,1,1 が答えになるし,もしもベクトルなら m,1 または 1,m が答えになります。
その他 この節では,基本的にベクトルを仮定しましたが,それは scilabでは 1行×m列 の行列として扱われていました。
なお,代入するときに要素をセミコロンで区切ることにより (例 a=[1;2;3;4;5]), n行×1列の配列も作ることができます。

操作の例:
CD に刻まれる音を題材にしてグラフ表示をしてみよう。
ここでは,1 ms の間の現象を確認したいと考える。 ここで,取り扱うべき変数 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);

これをプロットしてみる。

>> plot2d(t,f1);

カーソルを使えば元の行に直ぐに戻れる。 2kHz,10kHz,50kHz など幾つかの関数をプロットしてみよう。

以上は特に CD に限らない話であった。
ところで CD には,44.1 kHz という有限の時間でサンプリングされた情報が刻まれている。 (ここではグラフの書きやすさからサンプリング周波数を 40 kHz として話を進める。) すなわち,時間は飛び飛びになるのである。
fclk = 40 kHz ということは,
Tclk = 25 μs である。

先ほどの例で言えば,t を 250 個ごとに使うことに対応する。

それでは,サンプリングした情報をグラフ化してみよう。

>> plot2d(t(1:250:10001),f1(1:250:10001),style=-1);

注意:plot2d 関数の 3 番目の引数の style=-1 は,値をプロットすることを意味する
以上の操作のうち,打ち込む行を整理。ただし,f1 はどれか一つを選ぶ:
--> 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);
-->clf()
--> plot2d(t,f1);
--> plot2d(t(1:250:10001),f1(1:250:10001),style=-1);

あるいは次の行を実行しておくと,
>> tt = t(1:250:10001);
あとは,次の 1 行で全てをまかなうことができる。
>> f=1000; clf(); plot2d(t,sin(2*%pi*f*t)); plot2d(tt,sin(2*%pi*f*tt),style=-4)


演算c --- 文字列

文字列演算は,コンソール出力を行う際に必要な処理です。
代表的には以下のようなコマンドや,使い方を覚えておけば大丈夫。
文字列と入出力 ---きれいな出力を得るために---
<
概要 Scilab では,文字列は string型としてとして取り扱います。
代入操作例--> a='television' (文字列はシングルクォートで囲む)
string型操作例--> whos (上記の代入文の後なら,a は string 型)
(通常の数値は,constant 型)
操作例--> length(a) (文字数が得られる)
連結操作例--> b='This is a ' + a + '!!!' (加算記号で連結できる)
抜き出し操作例--> c=part(a,3:5) (3〜5文字目。上の例なら 'lev')


数値<==文字操作例--> x = 5 + evstr('123')
文字<==数値操作例--> msg = 'The answer is ' + string(x)


数値 x の表示操作例1--> x (次の行から3行かけて 'x = 128.' と表示
操作例2--> disp(x) (次の行から1行で '128.' と表示
操作例3--> disp(x,2x) (次の行から2行で '256. 128.' と表示
もしも一度に表示したい変数の数が,5個もあるならば, 操作例1 のように表示してると,数値以外が多くて見づらくなる。 こんなとき,disp という関数を使うと, 表示は変数の内容に限られるので,余分がなくて見やすくなる。 更に,上で示した文字列操作と絡めると,好みの出力が得られる。 なお,disp の引数は,文字列・数値・配列のいずれでも可能
応用例操作例--> disp( 'Data = ' + string(x) + ', Its squear =' + string(x * x) )



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

当然,2次元配列を扱うのも MATLAB の得意とするところです。
要素の取り扱い方法など,MATLAB 固有の簡潔な表記法があります。
代表的には以下のようなコマンドや,使い方を覚えておけば大丈夫。
演算---2次元配列,配列のグラフ化 (MATLABに違いがあれば)
概要 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 次元グラフ 
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 );
clf();
plot3d(x,y,z);
視点 Scilab を調査中

これ以降,配列専用の関数-2 (2次元配列)
zeros zeros(4,4) とは 4行4列 の,zeros(3,5) とは 3行5列の,全ての要素が 零 の配列
(MATLAB と 微妙な違いあり。焼き直しの際は注意が必要)
ones zeros と同様だが,全ての要素が 1 の配列
rand zeros と同様だが,それぞれの要素が 0 から 1 の間の乱数から成る配列
length(a) 配列 a について,要素の数,すなわち行の数×列の数 が返ってくる
(MATLAB と 微妙な違いあり。焼き直しの際は注意が必要)


演算e --- マトリクス

Scilab は,MATLAB がそうであるように,マトリクス演算を得意としています。
代表的には以下のようなコマンドや,使い方を覚えておけば大丈夫。
演算---マトリクス
概要 配列の演算はたいてい使える
行列の設定 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 の行と列を入れ替えたもの


プログラミング(*.sce-ファイル)

Scilab では,*.sce という拡張子を持つテキストファイルは, プログラムを格納するファイルとして使われます。関数は *.sci ファイルに格納します。
(MATLAB では .m ファイルです)
代表的には以下のようなコマンドや,使い方を覚えておけば大丈夫。
プログラミング(sce-ファイル)
プログラム実行の概要(Scilabのプログラム環境)
パス表示操作-->pwd (カレントディレクトリのパスを見る)
別の操作:メニューから「file」>「Get Current Directory」としても可能
パス変更操作-->chdir path (カレントディレクトリを,path に設定する)
別の操作:メニューから「file」>「Change Directory」として操作しても可能
ファイルのリスト操作-->dir (カレントディレクトリのファイルリストを,コンパクトに見る)
操作-->unix_x("dir") (カレントディレクトリのファイルリストを,詳細に見る)
エディタ起動 操作:メニューから「Editor」と操作する。
エディタ起動後に,操作すべきファイルを指定したり,セーブする
なお,Scilab では *.sce という拡張子が基本だが,拡張子をわざわざ .sce としてつけること。
特に拡張子を指定しないと,拡張のないファイルが出来てしまう。
sceファイルの編集 操作:メニューから「File」>「Open...」と操作し,ファイルを指定する。
別の操作:メニューから「Editor」と操作し,エディタの中からファイルを開く
実行操作:メニューから「File」>「Exec...」と操作し,ファイルを指定する。 別の操作-->exec mytest.sce (カレントディレクトリ内の mytest.sce というファイルに書かれたプログラムを実行させる)

sce-ファイルの作り方(Scilab プログラムの書き方)
概要 普通の演算はそのまま書いておけば大丈夫である
コメント 行のうち // (2つのスラッシュ)以降はコメントとみなされる
一時停止 pause; キーを押すまで一時停止
-1-> というプロンプトが現れて停止する。
継続には return を,終了には abort を打ち込む。
入出力 a=b+c; 値をあえて表示する必要が無い時
a=b+c デバッグのためにいちいち値を知りたいとき
disp(a) 変数名の表示は不用の時
disp( 'a =' + string(x) + ', Its squear =' + string(x * x) ) 細かい表示付き
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 行目参照)と, ファイル名は別で構わない。 続いて,
getf('myfunc2.sci');
というコマンドを実行すると,ファイル myfunc2.sci が読み込まれ,関数 myfunc1 が使えるようになる。 例えば y = myfunc1(5) によりこの関数が呼び出され, y は 120 になる。
<myfunc2.sci>
(Download)
function y=myfunc1(x)
s=1;
for n=2:x
 s=s*n;
end
y=s;
endfunction	


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

Scilab は複素計算もグラフ表示も簡単にできるため, 交流回路への応用も簡単です。
実行例(他にも良いやり方はいろいろあります)
応用---交流回路の電圧波形を計算(グラフの書き方)
系: R-L 直列回路
抵抗コイル入力電圧
定数値 200 Ohm1 H141 V

「50 Hzでの各部の波形をグラフ化するプログラム」

// s111
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);
clf();
plot2d(t,real(Vin),style=1)
plot2d(t,real(VR),style=2)
plot2d(t,real(VL),style=3)

「各部の電圧の周波数依存性を求めるプログラム」

// s112
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) ;
clf();
plot2d(f,abs(VR),logflag='ll',style=1)
(plot2d の引数 'll'(小文字の L が二つ)は, 縦軸も横軸も Log であることを表す。 'ln' は横軸だけが Log であり,'nl' は縦軸だけが Log である。 'nn' はデフォルト設定であり,両軸とも通常目盛りを意味する。)
本来は,VR と VL の2つのグラフを同時に書かせたいが, グラフの座標が変わってしまうので,ここでは一つずつ書かせる。 → VL も書かせてみよう


// s113
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) ;
clf();
plot2d(f,abs(VR),logflag='ll')
注意:ある数 (例えば f ) の逆数を求める際は,空白を入れながら書こう。 (例えば 1 ./ f のように。) 空白無しに (1./f) と書いたならば, (1. / f) = (1.0 / f) = マトリクスの割り算 と判断されるのか,又は, (1 ./ f) = 要素ごとの割り算 と判断されるか分らない。
注意: V はスカラーであるため(配列ではないので), 続く掛け算は * でも良いが, .* でもエラーにならないので, むしろ打ち間違いを防ぐために積極的に .* を使うことを勧める。


// s114
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 ;
clf();
plot2d(f,abs(VR2),logflag='ll')
注意: 並列は,a. * b. / (a + b) である。


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

MATLAB を使えば,電卓感覚で実験のデータ整理ができます。
実行例(他にも良いやり方はいろいろあります)
応用---実験データの整理(データと,回帰線)
例として,次の実験データについて回帰線を求める。
x(cm) 01245
y(cm) 0.10.92.89.015.2


// s121
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;
clf();
plot2d(x,y,style=-1);
plot2d(xx,yy,style=2);

(p1b1 で reglin(x,y) は1次式近似を求めます。 plot2d では style が負のときはプロットです。 )


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

MATLAB は信号処理も簡単に行なえます。
実行例(他にも良いやり方はいろいろあります)
応用---Waveファイルと信号処理(FFT,フィルタ)
もしもどこからでもいいから,*.wav というファイルを探し, ファイル名を sound.wav とし,カレントディレクトリに入れておくと, そのファイルに対してFFT処理を適用できます。


// s131 - FFT

Fs=44100
t=0:1/Fs:5;
y=0.05*sin(2*%pi*440*t)./(0.1+(t-0.3).^2) ...
 + 0.05*sin(2*%pi*770*t)./(0.1+(t-1.1).^2) ;

//[y,Fs,bits]=wavread('sound.wav');

playsnd(y,Fs,bits);

m=8192;
d=0;
clf();
for n=1:7

yy=y(1,1+(n-1)*m:n*m);
//wavsave('sound2.wav',yy);
playsnd(yy,Fs,bits);

scf(d);
// d=d+1;

subplot(2,1,1);
plot2d((1+(n-1)*m:n*m)/Fs, yy,style=n);

zz=fft(yy);

subplot(2,1,2);
plot2d((1:m/16)*Fs/m+60*n, ...
   abs(zz(1:m/16))-200*n,style=n);

end

//end

[y,Fs,bits]=wavread('sound.wav');
wavファイルを読み込みます。
y にはサンプルした音の各時間ごとの強度が入ります。 Fs には1秒あたりのサンプル数が, Bits には信号のビット数が入ります。 また信号がモノラルかステレオかにより y は1行n列か,または2行n列になります。
playsnd(yy,Fs,Bits);
yy を再生します。
y=0.05*sin ...
y を設定します。
savewave('sound2.wav',yy);
yy を保存します。
yy=y(1,1*Fs:1*Fs+8191);
y のうち,2秒目から 8192 点を抜き出します。 もしもCDと同じ44.1kSpsだったならば, 8192 点とは 約0.18秒です。
8192は,2の13乗であり,FFTをするのに最適です。
zz=fft(yy);
zz には yy のFFT結果が入ります。 yy が 8192 = 2^13 点 なので,高速にフーリエ変換出来ます。 zz もやはり 8192点 の配列になります。 サンプリング周期 Fs と,8192点であることから, abs(zz(1:4096))というデータは,周波数 Fs/2 [Hz] の 周波数成分の強度を表わします。
plot((1:4096)*Fs/8192, abs(zz(1:4096));
zz は複素数ですから,そのままではグラフ化出来ないので, 絶対値を取ってから表示します。



// s132 - Filter

Fs=44100;
bits=16;
t=0:1/Fs:5;
yin=0.3*sin(2*%pi*100*t)+0.3*sin(2*%pi*2000*t);

//[yin,Fs,bits]=wavread('sound.wav');

playsnd(yin,Fs,bits);

filt=[150 150]*2/Fs;
sys=iir(4,'lp','butt',filt,[0 0])

yout=flts(yin,sys);

playsnd(yout,Fs,bits);

n=1:2000;
plot2d(t(n),yin(n))
plot2d(t(n),yout(n),style=2)

[y,Fs,bits]=wavread('sound.wav');
y には音データが入ります。

filt=[150 150]*2/Fs;
sys=iir(4,'lp','butt',filt,[0 0])
フィルタを設定します。
4次です。
LowPassです。(他に'hp','bp','sb')
カットオフは 150Hz です。
[0 0]は'butt'には無意味です。

yout=filt(yin,sys);
先に設計した定数 sys に従って,
フィルタをyinに作用させ,
出力youtを得ます。
このプログラムは,次のウエブを参考にいたしました。
http://www.neurotraces.com/scilab/scilab2/node1.html
http://www.neurotraces.com/scilab/scilab2/node46.html


応用1d --- 数値積分

数値積分も MATLAB は色々な関数を取り揃えてあり,簡単に行なえます。
数値積分は,微分方程式の最も特殊な場合と考えることも出来ます。
実行例(他にも良いやり方はいろいろあります)
応用---数値積分
まずは,以下の関数を作り,myfunc3.m という名前で保存(Save)しておく。 これでファンクションが定義できたことになり,以降たとえば,
y=myfunc3(3.14159/2)
とでも打ち込むと,y にはほぼ 1 (= sin(90度) )が設定される。
実際のところ積分の計算には必ずしもfunctionが必須というわけではないが, こうすることによって(1)プログラムが見やすくなる
(2)function を使うと大きな利点もある(最後の例を参照のこと)ので,お勧めである(というか,functionの利用を前提として説明する)。
<myfunc3.sci>
(Download)
function y=myfunc3(x)
y=sin(x);
endfunction

さて,まずは自分でプログラムを作り, 区間 0〜%pi にわたって積分しよう。 採用するアルゴリズムは,一番簡単な台形法とする。
T=h*(y1 + 2*y2 + 2*y3 + yn)/2

// s141
// clear myfunc3
getf('myfunc3.sci');
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))
(p1d1 の clear 文は,メモリ内に残る myfunc3 を消す。
これは,myfunc3 を別の関数に変えたときに, メモリ内に残る myfunc3 を消すためである。 なお,clear 文の前の // は最初から消してしまってかまわない。 というのも, myfunc3 の読み込み前に clear myfunc3 を行ってもエラーにならないからである。)

上記の台形法と同じアルゴリズムで計算する関数が用意されており, 次のように行っても同じ結果が出る。
% p142 - MATLAB による積分の関数(アルゴリズムは台形法)を利用
n=201;
x=linspace(0,%pi,n);
y=myfunc3(x);
s=trapz(x,y)
% end of program
(Scilab については調査中)

functionの形で積分関数が設定されている場合, 組み込み関数を使うと,積分はたったの一命令で実行できる。
組み込み関数による積分
// s143 
s=intc(0,%pi,myfunc3)
注意: MATLAB と Scilab を比較すると, 積分を行う関数は,名前が違うだけでなく, パラメタを与える順番も,パラメタの与え方 ('...クォーテーションマーク で囲むかどうか)まで異なるので, プログラムを移転する場合は注意が必要である。


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

微分方程式を考える上で,2次の微分方程式は基本中の基本で, 応用範囲の広いものです。 それは,変位と速度から加速度を導く系になぞらえることが出来ます。 また,RLC回路の応答にもなぞらえることが出来ます。
そういうわけで,ここでは,2次の微分方程式の解法を例に取り, MATLAB で解く方法を考えたいと思います。
取り扱う系は,振り子の振動を表わす,
x'' = - a x' - sin(x)
とします。
実行例(他にも良いやり方はいろいろあります)
応用---微分方程式
ここでは,次の式の解を計算によって求めて,振り子の振動を解析するものとする。
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
となる。この式は,左辺は 1 次の微分だけであり, 右辺は通常の項だけから成る。
これは以下のようにプログラムにできる。 なお,今回は関数内で t を使っていないが, 一般に微分方程式は時間の関数でもあることから, 引数の中には時間 t も入れておく。
<myfunc4.sci> (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

さて,まずはこの微分方程式を解くにあたって, まず 変位-速度 平面上での解軌跡を考える。 以下のプログラムを打ち込んで実行してみよう。 この図表から解のふるまいがかなり予想できる。
変位-速度平面上での解軌跡
横軸を変位 x ... x2,
縦軸を速度 v ... x1 とする
注釈

// s151 
// x...x2, v...x1

clear myfunc4
getf('myfunc4.sci');

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

clf(); champ(x, v, dx(:,:,2), dx(:,:,1))
最初に定義した  x ベクトルと y ベクトル:
x=-9:0.5:9;
v=-5:0.5:5;
は,グラフの横軸と縦軸の変域を指定するものである。
グラフを拡大・縮小するには範囲やステップを調整すればよい。

myfunc4 を呼び出すときのパラメタは [v(nv) x(nx)] である。
このことから,
・パラメタベクトルの最初の要素は v に関連し, 
・次の要素は x に関連するといえる。
これより,myfunc4 は次のように解釈できる。

function xdot=myfunc4(t,x)
 xdot=zeros(2,1); 
 V=x(1); X=x(2); 
 a=0.5;
 V_dot =-a*V-sin(X);
 X_dot = V; 
 xdot(1)=V_dot; xdot(2)=X_dot;
endfunction
変数名を変更すると,印象が変わります。
横軸は変位 x ... x2,
縦軸は速度 v ... x1.

// s151 - ver2 
// x...x2, v...x1
clear myfunc4
getf('myfunc4.sci');
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
clf(); champ(x1, x2, dx(:,:,2), dx(:,:,1))
MATLAB と Scilab を比べると, 同じ機能の関数同士でもパラメタの与え方が異なることがある。 ここに示した例では,出力した図は同じであるにもかかわらず, dx というマトリクスのサイズは MATLAB では 21 * 37 * 2 であるのに対し, Scilab では 37 * 21 * 2 である。

それでは,この微分方程式を, 組み込み関数を使わずに解いてみる。 今回は横軸を時刻としてグラフを書く。
振り子の振動

// s152 
// x...x2, v...x1

clear myfunc4
getf('myfunc4.sci');

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

clf();
plot2d(tt,xx(:,1),1);
plot2d(tt,xx(:,2),2);

微分方程式を解くための関数 ode を使うと, このような微分方程式を簡単に解くことができる。 早速"help ode"と打ってみよう。
なお,ここでは,グラフを2回書く。 一枚目は p1e2 と同じである。 二枚目は横軸を変位,縦軸を速度としたもので, p1e1 に対応するものである。
一枚目のグラフを書いたら pause で停止するので, return(c/r) を打って,二枚目のグラフに進む。
振り子の振動を時間軸で解く

時間領域
: 0〜20 を 0.1刻み(秒)

ODEは微分方程式を解く関数
◆[4;0] とは,初期条件
 x(1)=V=4;
 x(2)=X=0;
◆0 とは,初期時間
◆tt は解を求める時刻ベクトル。
必ずしも 0 を含む必要はない。
◆myfunc4 とは解く方程式を
プログラムした関数

// s153
// x...x2, v...x1
clear myfunc4
getf('myfunc4.sci');
tt=0:0.1:20;
xx=ode([4 ; 0], 0, tt, myfunc4);
clf();
plot2d(tt,xx(1,:),1);
plot2d(tt,xx(2,:),2);
pause;
clf();
plot2d(xx(2,:),xx(1,:));


せっかく解いたので, 最初に示した平面上に解軌跡を表示させてみる。 Scilab では特に画面を消去しない限り同じグラフ画面上に重ねがき出来るので, 別々に求めた解を,同じグラフ上に表示する。
変位-速度平面上での
流れと解軌跡

x ... x2,
v ... x1.

プログラム s151
プログラム s153
を合体させたもの。

// s154

clear myfunc4
getf('myfunc4.sci');
clf();

// from s151
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 s153
tt=0:0.1:20;
xx=ode([4 ; 0], 0, tt, myfunc4);
plot2d(xx(2,:),xx(1,:),2);


以上のように,この節では微分方程式の解き方について, ユーザーがすべて解くやり方と,組み込み関数(ode)を使った場合の 2 通りの解き方を示した。
自分で方程式を解く方法は,微分方程式の考え方の習得に有効であった。
一方,組み込み関数を使う方法は,極めて少ないプログラミングでも, 実用的な問題解決に有効であることを示した。
ここで示さなかった微分方程式に対しても同様に取り組むことが出来るので参考にして欲しい。


応用---微分方程式-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, p1e4 と同じであり, 変更点は色違い(赤)で示した。 ただし,コメント行については変更があっても特に示してない。

<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
% p1e5 - 角度-角速度上の解軌跡と,時間-角度,角速度のグラフ

%       横軸 は 角度 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 (符号) 関数で与えられるため,方形波である。 従って,2通りの力の掛かり方しかないため,ベクトルの図は必要十分の2枚書くことにした。
2つのベクトルの方向が同一の場合は, 後から書いた矢印によって,最初に書いた矢印の一部が消されているので, しまうので注意しよう。

ここで,上記とは異なる非常に長い周期の外力を加えた場合の解を示す。 二通りの外力に応じた二通り場所に収束しているのがわかる。
myfunc5, p1e5 からの変更点は次の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, p1e5 と同じであり, 変更点は色違い(赤)で示した。 ただし,コメント行については変更があっても特に示してない。
<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.sci> (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)
// s161
clear myfunc7
getf('myfunc7.sci');
tt=0:1.e-4:10.e-3;
yy=ode([0], 0, tt, myfunc7);
vc=yy(1,:);
clf();
plot2d(tt,vc,1);

RC 回路の時間応答 (Rの電圧)

同じ回路を用いて, 抵抗にかかる電圧 Vr を出力とした場合についてグラフ化しよう。 式は,
 
Vr = τ s
1 + τ s
Vin = { 1 - 1
1 + τ s
} Vin
で表わすことができる。ただし,τ=1/(CR) である。
二つ目の等号から,Vr ← Vin - Vc であることがわかる。
従って,p1f1 の続きとして, 次の命令を続ければ,Vr をグラフ化することができる。

% prog1f1 の続き
vr=1-vc;
plot(tt,vc,tt,vr)
// s161 の続き(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の過渡応答を解釈するには, 2つの時定数を考えればよい。 一つは,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 の特性と言える。

先のp1f1 にて VC2 は既にグラフ化した。
また,先立つ p1f1 の「経過2」の計算から, s * VC2 は y2 即ち,yy(:,2) として求められる。

従って,p1f2 の続きとして, 次の命令を続ければ,VR2 をグラフ化することができる。
% prog1f2 の続き
vr2=100.e3 * 1.e-6 * yy(:,2);
plot(tt,vc2,tt,vr2)
// s162 の続き(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 乗の項を引くことである。
従って,p1f2 の続きとして, 次の命令を続ければ,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
// s162 の続き<開発中>

※ 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
で表わすことができる。
先立つ p1f3 の「経過」の計算から, s * VC は y2 即ち,yy(:,2) として求められる。
% prog1f3 の続き
hold on
vr=1.e-9 * 1000 * yy(:,2);
plot(tt,vr,'B'); % B---Blue
// s163 の続き(Scilab)<開発中>

続いて,出力を同じ回路のコイル L として計算する。 この項を求めるには,1 (=Vin) から,抵抗やコンデンサの電圧を引くことである。
% prog1f3 の続き
vL = 1 - vc - vr;
plot(tt,vL,'R'); % R---Red
hold off
// s163 の続き(Scilab)<開発中>

※ 確認しよう : R の値を変えると,減衰率が変化する。 CまたはLを変えると振動の周期が変化する。

※ Vin がランプ波(1次式で時間と共に増加する波形)だった場合, は,本来 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 を使うのであって, 各項の大きさについては特に物理的な根拠はありません。

方程式:n型半導体内の少数キャリア(正孔)の運動を表す式

(ここで, ' とは,x による微分を表す)

dp/dt = D * p'' - u * E * p' - u * p * E' + G - (p - p0)/T

ただし,
p は正孔濃度
p0 は熱平衡時の正孔濃度
D は拡散定数
u は移動度
E は電界
G は生成項
T は再結合時間

まずは,このプログラムが動くように設定し, 一度走らせて見よう。 解説については,その後に行う。


半導体中の正孔の振舞

 // s171
 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;
 
 clf();
 plot2d(x,p);
 plot2d(x,4.5+ee,style=3);
 ee=0.02*ee;
 plot2d(x,6+g,style=4);
 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,style=2,rect=[0 0 1400 8]);
 disp('--end--');
 // end of program



それでは,この微分方程式を微妙に変えてみよう。
  • 注意:この計算の中で,ドリフト電流を計算する際には, 拡散項は省くことができない。
  • // 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 行は, それぞれ生成項と,生滅項の計算部である。 これらの行をコメントアウトすれば, その物理現象のシミュレーションを省くことができる。