• プログラムの実行方法は Scilab超入門 で独習できます.
  • プログラム例は,SciNotes にコピー&ペーストして使ってください.

もどる

12章と15章

//
// "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

14章

//
// "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();

もどる