- プログラムの実行方法は Scilab超入門 で独習できます．
- プログラム例は，SciNotes にコピー＆ペーストして使ってください．
テキスト12章のプログラム例
//
// "pend.sce" //から右は処理されない．「コメントアウト」というメモ用の文法
//
clear; clf(); //変数リセット; グラフリセット;
//運動方程式を解く
M=2/3; m=1/3; l=1;
function dx = eom(t,x) //運動方程式の定義
ft = 0.0; //← (1)ここに制御力！
A = [M+m, m*l*cos(x(3)); cos(x(3)), l];
b = [m*l*(x(4)**2)*sin(x(3))+ft; 9.8*sin(x(3))];
h = A\b; //(Aの逆行列)*bと同じ．inv(A)*bと書くより速い
dx(1) = x(2); dx(2) = h(1);
dx(3) = x(4); dx(4) = h(2);
endfunction
x0=[0; 0; 0.5; 0]; // ← (2)ここに初期値！
n=200; tt=linspace(0,25,n); //時間を表す等差数列
xx=ode( x0, 0, tt, eom ); //運動方程式を解く
//アニメーションする．
function draw_mech(x)
g=gca(); g.isoview="on"; //座標軸の取得; 縦横比1;
g.data_bounds=[-4,-1.5;4,1.5]; //座標軸の設定
xM = [x(1);0];
xm = xM + l*[sin(x(3));cos(x(3))];
plot([xM(1),xm(1)],[xM(2),xm(2)],'r-'); //振り子
p=gce();p.children.thickness=3; //直前の描画の線の太さ
xlabel("xx21xx"); //学籍番号
xgrid(2); //グリッドon
endfunction
draw_mech(x0); drawnow; //初期描画; 画面更新;
sleep(2000); //2秒待ち
realtimeinit(0.1); //アニメーションの時間刻み
for i=1:n //コマ送り
realtime(i); //リアルタイム更新設定
drawlater(); clf(); //描画延期; 画面消去
x=xx(:,i); //状態ベクトル取得
draw_mech(x); drawnow; //機構描画; 画面更新
end
旧テキストのプログラム例
//
// "vib-tr.sce"
//
clear; clf();
function dx = model(t,x)
om=1.6;
dx(1) = x(2);
dx(2) = -0.2*x(2) - x(1) + cos(om*t);
endfunction
x0 = [0; 0.3]; tt = linspace(0, 100, 800);
xx = ode(x0, 0, tt, model);
g=gca(); g.data_bounds=[0,-1.5;100,1.5]; //座標軸の設定
plot(tt,xx(1,:),"-");
xtitle("x(t)=a(t)+b(t)"); xgrid();
xlabel("t"); ylabel("x(t)");
//
// "vib-res.sce"
//
clear; clf();
function dx = model(t,x)
dx(1) = x(2);
dx(2) = -0.2*x(2) - x(1) + cos(om*t);
endfunction
om1 = linspace(0.2,1.6,15);
x0 = [0; 0.1]; tt = linspace(0, 100, 300);
realtimeinit(0.1); //アニメーションの時間刻み
for i = 1:15
realtime(i);
drawlater(); //描画延期
om = om1(i); f = cos(om*tt);
xx = ode(x0, 0, tt, model);
clf(); subplot(2,1,1);
plot(tt,f,"-",tt,xx(1,:),"-");
xlabel("t"); ylabel("x(t)");
xtitle(sprintf("Waveform (om = %.1f)",om));
g=gca(); g.data_bounds=[0,-5;100,5]; xgrid(); //座標軸の設定
drawnow; //画面更新
xxmax(i)=max(xx(1, 250:300)); //xx(250,1)〜xx(300,1)の最大値
end
subplot(2,1,2); plot(om1, xxmax, "o-");
xlabel("om"); ylabel("max x(t)");
xtitle("Response Curve");
save("vib_res.dat","om1","xxmax"); // vib-res2.sce で使うデータの保存
//
// "vib-res2.sce"
//
clear; clf();
load("vib_res.dat","om1","xxmax"); // vib-res.sce で保存したデータ
m=1; c=0.2; k=1; P=1;
zeta = c/( 2*sqrt(m*k) ); // 表3.1
omn = sqrt(k/m); // 表3.1
function y = K(Om)
global zeta;
y=1/sqrt((1-Om^2)^2 +(2*zeta*Om)^2);
endfunction
om2 = linspace(0.2,1.6,100);
A = P/(m*omn^2);
for i=1:100
R(i) = A*K(om2(i)/omn);
end
plot(om1, xxmax,"o", om2, R,"-" );
xlabel("om"); ylabel("Amplitude");
xtitle("Response Curve (o max(x); - R)"); xgrid();
g=gca(); g.data_bounds=[0.2,0;1.6,6]; xgrid(); //座標軸の設定
//
// "vib-bode.sce"
//
clear; clf();
m=1; c=0.2; k=1; P=1;
zeta = c/( 2*sqrt(m*k) ); //表3.1
omn = sqrt(k/m); //表3.1
function y = K(Om)
global zeta;
y=1/sqrt((1-Om^2)^2 +(2*zeta*Om)^2);
endfunction
function y = phase_diff(Om)
global zeta;
y=-atan((2*zeta*Om),(1-Om^2)); //atan2(y,x)
endfunction
om = linspace(0.2,1.6,100);
A = P/(m*omn^2);
for i=1:100
R(i) = A*K(om(i)/omn);
phi(i) = phase_diff(om(i)/omn);
end
subplot(2,1,1); plot(om, R);
xlabel("om"); ylabel("Amplitude");
title("Amplitude Ratio"); xgrid();
subplot(2,1,2); plot(om, phi);
xlabel("om"); ylabel("phi");
title("Phase Difference"); xgrid();