[1] 行列データの入力
1)サイズ2×2の行列
を定義するために、次を入力してみます。
>> A=[1 2;3 4]
ここで、行列データは[ ]で囲み、要素の区切りはスペースまたはカンマ、行替えはセミコロン;を用いることに注意してください。次のように最後に;を置くと、出力が抑制されます。
>> A=[1 2;3 4];
すでに入力した行列データへのアクセスは、次のように行います。
>> A
また、特定の要素へのアクセスや修正は、次のように行います。
>> A(2,1)
>> A(2,1)=-5
2)2次の単位行列
は、次により入力できます。
>> B=eye(2,2)
便利なのは、次のように異なる行数と列数が許されることです。
>> B=eye(2,3)
これは、次のサイズ2×3の行列
の定義を意味します。
3)サイズ3×2の零行列
は、次により定義できます。
>> C=zeros(3,2)
4)3次の対角行列
は、次により定義できます。
>> D=diag([1 2 3])
すなわち、ベクトルにdiagを作用させているわけですが、行列にdiagを作用させ、対角成分を抽出することもできます。
>> d=diag(D)
5)行列A,B,C,Dを小行列としてもつブロック行列
を定義するためには、次のように入力します。
>> E=[A B;C D]
6)サイズ2×2の複素行列
も、次により定義できます。
>> F=[1+i, 2+3*i;3+2*i, 4+4*i]
ここで、iは虚数単位として定義されています。もちろん再定義可能です。
7)MATLABでは、空行列を定義できます。
>> G=[]
8)すべての要素が1である行列
は、次により定義できます。
>> H=ones(2,2)
[2] 行列データの管理
1)現時点でどのような行列を定義したかを調べるためには、次のように入力します。
>> who
行列のザイズを合わせて表示するには、次のように入力します。
>> whos
また、たとえば入力した行列データdを消去するためには、次のように入力します。
>>clear d
すべてのデータを消去するときは、次のように入力します。
>>clear
2)行列データEのサイズを調べるには、次のように入力します。
>> [m,n]=size(E)
3)浮動小数点表示を行うためには、次のように入力します。
>> format short e
これを固定小数点表示に戻すためには、次のように入力します。
>> format
4)カレントディレクトリのチェックを行うためには、次のコマンドが使えます
>> pwd
>> ls
カレントディレクトリの変更のためには、次のコマンドを与えます。
>> cd dname
5)メモリ上のデータのセーブとロードのためには、次のコマンドを与えます。
>> save filename
>> load filename
[3] コロン記法
1)まず、次を入力してみてください。
>> h=1:3
これは、「form 1 to 3」を意味します。すなわちサイズ1×5の行列h=[1 2 3]が生成されます。コロン「:」はtoの意味で使われ、colon notationと呼ばれています。
2)それでは、次はどうでしょう。
>> h=0:0.5:3
これは、「form 1 to 3 by 0.5」を意味します。すなわち、2つのコロンで挟まれた0.5は増分を意味し、サイズ1×5の行列h=[1 1.5 2 2.5 3]が生成されます。ちなみに、h=1:3はh=1:1:3と同じと考えておきましょう。
3)それでは、次はどうでしょう。
>> h=0:-2
これは、h=1:1:-2の意味になりますが、1から-2まで、増分1では到達できないので、空行列となります。
[4] 小行列の抽出
1)コロン記法の利点は,行列の添え字部分に適用できることにあります。たとえば,行列Eの(2,3)要素には
>> E(2,3)
でアクセスできるのですが,この2と3は,スカラーである必要はなく,ベクトルが使えます。たとえば,次を確かめてください。
>> E([1 3],: )
これは第1行と弟2行からなるサイズ2×5の小行列が取り出されることに注意してください。列を指定する「:」はすべての列を意味しています。
2)次はどうでしょう。
>> E(:,2:4 )
ここで、2:4は[2 3 4]というベクトルでしたよね。第2~4列からなるサイズ5×3の小行列が取り出されることに注意してください。行を指定する「:」はすべての行を意味しています。
3)それでは,次はどうでしょう。
>> E([1 3],2:4 )
[5] 行列の分割
1)いろいろな応用では、ある行列を分割する必要が出てきます。これは小行列の抽出の観点から、たとえば次のようにできます。
>> [m,n]=size(E); i=2; j=2;
>> A=E(1:i,1:j);
>> B=E(1:i,j+1:n);
>> C=E(i+1:m,1:j);
>> D=E(i+1:m,j+1:n);
ここで、i=0のときは1:i=[]、i=mのときはi+1:m=[]などに注意してください。
[6] 行列演算
いま、行列Aが次のように定義されていることを確認してください。
また、行列Bを新しく定義します。
>> B=eye(A)
1)行列A,Bの和、差、乗算は次のように計算できます。
>> X=A+B
>> X=A-B
>> X=A*B
2)行列Aが正則のとき、その逆行列は
>> X=inv(A)
で計算されます。しかし、これは逆行列の要素が欲しいとき以外は使わないようにしましょう。
多くの場合、逆行列は次の連立一次方程式
または
を解くときに用いますが、MATLABでは、次のように計算します。
>> X=A\B
>> X=B/A
それぞれ左からの割り算、右からの割り算と呼ばれます。
これは数値的に安定な数値計算アルゴリズムを得るためは直交行列だけの演算ですませることが必要だからです。
3)行列Aのべき乗と転置共役は次のように計算でされます。
>> X=A^2
>> X=A’
4) それでは、次の演算結果と比較してみましょう。
>> X=A.+B
>> X=A.-B
>> X=A.*B
>> X=A.\B
>> X=B./A
>> X=A.^2
>> X=A.’
すなわち、演算子の前に「ドット.」をおくと、要素ごとの演算やただの転置になります。
5)次に、組み込み関数としては、次が使えます。
sin,cos,tan,asin,acos,atan,sinh,cosh,tanh,asinh,acosh,atanh,sqrt,exp,log,log10,abs,real,imag
6)行列関数としては次が準備されています。
>> X=expm(A)
>> X=norm(A)
>> X=rank(A)
>> X=det(A)
ぞれぞれ、行列指数関数、ノルム、階数、行列式を計算します。
7)また、固有値問題を解くには次の関数を用います。
–>[V,R]=eig(A)
単に固有値だけを求めるのであれば
–>r=eig(A)
でOKです。
8)さらに、特異値問題を解くには次の関数を用います。
–>[U,S,V]=svd(A)
単に特異値だけを求めるのであれば
–>s=svd(A)
でOKです。
9)多項式との掛け算は、次のように行います。
–>r=convol([1 1],[1 2])
10) いま、次のベクトルを定義します。
–>x=[3 1 5 4 2]
このとき、次はxの要素を小さい順番に並べ替えます。
–>[y,i]=sort(x)
yがソーティングされた結果、iが元のどの位置にあったかを示すインデックスです。たとえば
–>x(i)
としてみるとyと同じ結果が出ます。
11)次は、ベクトルの要素和、最大値、最小値を求める関数です。
–>y=sum(x)
–>y=max(x)
–>y=min(x)
[7] 関数のグラフ
1)次の関数のグラフを描いてみます。
次のコマンドを与えてください。
>> x=-3:0.1:3;
>> y=1/sqrt(2*pi)*exp(-x.*x/2);
>> plot(x,y)
ここで,なぜx*xでなく,x.*xとしているか考えてみてください。またpiはとして定義されています。
2)描画範囲(-3<=x<=3, 0<=y<=0.4)を指定し,グリッドを描画するには,次のようなコマンドを追加します。
>> grid
>> axis([-3 3 0 0.5])
>> xtitle(’N(0,1)’)
3) 次の1次系の周波数応答を片対数グラフび描いてみましょう。
次のコマンドを与えてください。
>> w=logspace(-1,1);
>> ga=20*log10(sqrt(1+w.^2).^(-1));
>> semilogx(w,ga)
>> grid
[8] プログラミング
これまでは対話型入力を行ってきましたが,一連のコマンドをエディタを用いてあらかじめ編集しておいて一括処理する,すなわちMATLABでのプログラミングについて学びます。まず,新規スクリプトを示すアイコンをクリックして,エディタを起動してください。
1) 1から9までの奇数の和を求める次のプログラムを入力し,ファイル名ex1.mでセーブします。
//ex1.m
//—–
s=0;
for i=1:2:9,
s=s+i;
end
disp(s,’sum=’)
//—–
//eof
エディタの実行ボタンをクリックします。答えは25のはずです。
2) マシンイプシロン(machine epsilon)と呼ばれる特別な定数を求めてみましょう。次のプログラムをファイル名ex2.mでセーブします。
//ex2.m
//—–
s=1;
while 1+s>1,
s=s/2;
end
s=s*2;
disp(s,’eps=’),
//—–
//eof
エディタの実行ボタンをクリックします。このsは次のepsと等しくなります。
>> eps
次を実行してみてください。
>> 1+eps/2==1
すなわち,
マシンイプシロンeps=a^{-k}は,1+a^{-k-1}>1が成り立たなくなる最小のkを用いて定義されます。浮動小数点による数値計算では,1+eps=1となる正数が存在することに注意してください。
3)次は,if文を使った例です。次のプログラムをファイル名ex3.mでセーブします。
//ex3.m
//—–
x=input(’x= ‘);
if imag(x)~=0,
disp(’complex number’)
elseif x>0,
disp(’positive real number’)
elseif x<0,
disp(’negative real number’)
else
disp(’zero’)
end
//—–
//eof
エディタの実行ボタンをクリックします。。
4)関数の例として,行列の階数を計算するものを示しましょう。次のプログラムをファイル名rank2.mでセーブします。
//rank2.m
//—–
function r=rank2(A)
s=svd(A);
tol=max(size(A))*s(1)*eps; r=sum(s>tol);
endfunction
//—–
//eof
このとき,次を入力してください。
>> r=rank([1 2;3 4])
5) 次の微分方程式の解を求めることを考えます。
次のプログラムをファイル名ex4.mでセーブします。
//ex4.m
//—–
function dx=f(t,x),dx=-x+sin(10*t), endfunction
x0=1; t0=0;
t=0:0.01:5;
x=ode(x0,t0,t,f);
plot(t,x),
mtlb_grid
//—–
//eof
エディタの実行ボタンをクリックします。