2020年3月22日日曜日

MATLABことはじめ 1

計算ツールはMathematica派だし,必要なものはPythonなりでコードを書いちゃえばいい派なので避けてきたのですが,制御をやる上でちょっと使ってみよう(使った上でPythonにもどるとか)と思い立ち始める次第です.

このエントリーは何かをまとめるというよりは,随時メモを乱投するスタンスで書いてみようと思います.

エントリーを区切ったり,追加したりは気分次第です

2020/05/10更新 微分方程式関連を追記

MATLAB

MATLAB = Matrix Laboratory の略とな

いろいろ楽をするためにできたものだと思うけど,過去資産を利用するための学習コストが見合うのか.いや,迷わず信じる.信じたものしか救われないのではないだろうか.私の信仰心がどこまで続くかを楽しんでもらえれば良いと思います.

各行の改行は;でおこなう. C言語とおなじ.

行列 

インデックスは0からではなく1から始まる. スペースで区切ってもカンマで区切ってもあり(oo).

行列の場合2行目以降は;で改行する.

{}はCell Array. 配列の配列を使用する.

a=[1 2 3]
b=[1,2,4]
c=[1,2;3,4]
%e=[a,b,c];
%disp(e); これはディメンジョンが合わないのでエラーになる
f={a,b,c};
disp(f);

配列要素へのアクセスは[]ではなく()を使う(oo).

% disp(c[1,1]);これはエラー
disp(c(1,1));

行列の演算.

x,y,zの座標を格納した連続データに対して指定した座標分平行移動したい.

次元があっていれば,各要素要素にたいして演算をしてくれる. 行と列を区別することに注意する


points(1,1) =  0;points(2,1) = 0;points(3,1) = 0;
points(1,2) =  5;points(2,2) = 0;points(3,2) = 5;
points(1,3) = 10;points(2,3) = 0;points(3,3) = 0.5;
points(1,4) = 15;points(2,4) = 0;points(3,4) = 5;
disp(points);

TranslationVector=[100;10;30];
%これをTranslationVector=[100,10,30];とすると次の演算でエラーになる
d=points+TranslationVector;
disp(d)
%やっていることは次の文と同じ. forはできるだけ使わないようにしたい(心意気)
for i=1:4
    C(1:3,i)=points(1:3,i)+TranslationVector;
end
disp(C);

微分方程式関連

1階の微分方程式
% 微分方程式をとく
% 時間tの関数, aは微分方程式の係数, dはt=0でのyの初期値
syms y(t) a d
% 微分方程式を定義
ode  = diff(y,t) == -a*y;
% 微分方程式を解く. 初期条件を指定しない場合
S=dsolve(ode)
% 微分方程式を解く. 初期条件を指定する場合
cond = y(0) == d;
dsolve(ode,cond)

% パラメータをあとから指定したい場合関数を定義して結果を代入
S2(t,a,d)=dsolve(ode,cond);
% a=2, d=1の例
S2(t,2,1)

%グラフ化してみる
% xを指定 初期値:ステップ:終値
t = 0:0.1:5;
%plot はx,y,線種を並べる
plot(t,S2(t,-1,1),t,S2(t,1,1),'--',t,S2(t,0.5,1),':')
% プロットの範囲指定
xlim([0 5])
ylim([0 2])
%凡例指定
legend('a=-1','a=1','a=0.5')
2階の微分方程式
%2階微分方程式を解く
% 時間tの関数, a1, a2は微分方程式の係数, dはt=0でのyの初期値
syms y(t) a1 a2 d
dy = diff(y,t,1);
% 微分方程式を定義
ode  = diff(y,t,2) +a1*dy + a2*y == d;
% 初期条件は低次からリストの形で与える
cond = [y(0)==1, dy(0)==0];
% 初期条件をつかわない場合
S3(t)=dsolve(ode)
% 初期条件を使って解く場合
S4(t,a1,a2,d)=dsolve(ode,cond)

%グラフ化してみる
% xを指定 初期値:ステップ:終値
t = 0:0.1:5;
%plot はx,y,線種を並べる
plot(t,S4(t,11,10,1),t,S4(t,2,10,1))
% プロットの範囲指定
xlim([0 5])
ylim([-1 1])

0 件のコメント: