Octave 3.0.5 (Linux) で動作確認したプログラム例です．

実行環境によっては，簡易アニメが滑らかに動かない場合があるそうです（plot の仕様変更？）．実行環境の違いにより，皆様の不具合を当方が再現できず，有効な対策をご提案できないケースが起こりそうですが，その点ご容赦下さい． m(_ _)m

Code 1 “pend.m”

# "pend.m"　　ここはコメント行．#～行末をOctaveは無視する
global m l c g;  #大域変数の宣言
m=1; l=1; c=1; g=9.8;
function dx=eqn(x,t)
 global m l c g; #関数の中から大域変数を参照する宣言
 dx(1) = x(2);
 dx(2) = -c/(m*l*l)*x(2)-g/l*sin(x(1));
endfunction
x0=[pi/5;0];             #初期値．pi は円周率
n=100;                   #時間のきざみ数
tt=linspace(0,10,n);     #時間0～10の等比数列
xx=lsode("eqn", x0, tt); #常微分方程式"eqn"を解く
for i=1:n  #簡易アニメーション
 theta=xx(i,1);          #第1変数が角度
 plot( [0,sin(theta)], [0,-cos(theta)], 'linewidth', 4); #線分のプロット
 axis([-2,2,-2,2],'square');    #座標軸の設定
 title(sprintf("xxxxx: %d",i)); #タイトルの設定
 grid on; drawnow();
 if ( i==1 ) sleep(1);   #初回は1秒待機
 else sleep(0.1); endif  #この数値でアニメーション速度を調整
endfor

Code 2 “vectorfield.m”

※Code 3, 4 の実行に必要．

# "vectorfield.m"
function vectorfield(func,nx,ny,rx,ry,scal)
 xx=linspace(rx(1),rx(2),nx);
 yy=linspace(ry(1),ry(2),ny);
 uu=zeros(ny,nx); vv=zeros(ny,nx);
 for ix=1:nx
  for iy=1:ny
   x = [xx(ix),yy(iy)];
   eval(sprintf("dx=%s([xx(ix),yy(iy)],0);",func));#dx = eqn(x,0);
   uu(iy,ix) = dx(1); vv(iy,ix) = dx(2);
  endfor
 endfor
 quiver(xx, yy, uu, vv, scal, "k-");
endfunction

Code 3 “pendvec.m”

※同じフォルダに，Code 2を書き込んだファイル vectorfield.m が必要．

# "pendvec.m"
global m l c g;
m=1; l=1; c=1; g=9.8;
function dx=eom(x,t)
 global m l c g;
 dx(1) = x(2);
 dx(2) = -c/(m*l*l)*x(2)-g/l*sin(x(1));
endfunction
nx=20; ny=15; rx=[-2,8]; ry=[-6,6];
hold off;
vectorfield("eom",nx,ny,rx,ry,0.5);
axis([rx(1),rx(2),ry(1),ry(2)]);
tt=linspace(0,10,200);
hold on;
xx=lsode("eom",[1.2;6],tt);
plot(xx(:,1),xx(:,2),"k-",'linewidth',5)
xx=lsode("eom",[1.4;6],tt);
plot(xx(:,1),xx(:,2),"k-",'linewidth',5)
plot([pi],[0],"ko",'markersize',10);
xlabel("x",'fontsize',32); ylabel("dx/dt",'fontsize',32); drawnow();

Code 4 “pendvecL.m”

※同じフォルダに，Code 2を書き込んだファイル vectorfield.m が必要．

# "pendvecL.m"
global m l c g;
m=1; l=1; c=1; g=9.8;
function dv=eom_v(v,t)
 global m l c g;
 dv(1) = v(2);
 dv(2) = -c/(m*l*l)*v(2)-g/l*v(1);
endfunction
function dw=eom_w(w,t)
 global m l c g;
 dw(1) = w(2);
 dw(2) = -c/(m*l*l)*w(2)+g/l*w(1);
endfunction
nx=20; ny=15;
clf; hold off;
rx=[-2,8];ry=[-6,6];
vectorfield("eom_v",nx,ny,rx,ry,0.5);
axis([rx(1),rx(2),ry(1),ry(2)]);
tt=linspace(0,10,200);
hold on;
xx=lsode("eom_v",[0.5;-4.8],tt); #[1.2;6]
plot([0.5],[-4.8],"ko",'markersize',10);
plot(xx(:,1),xx(:,2),"k-",'linewidth',5)
xlabel("v",'fontsize',42); ylabel("dv/dt",'fontsize',42); drawnow();
sleep(3);
clf; hold off;
rx=[-2-pi,8-pi];ry=[-6,6];
vectorfield("eom_w",nx,ny,rx,ry,0.5);
axis([rx(1),rx(2),ry(1),ry(2)]);
hold on;
plot([1.4-pi,1.6-pi],[6,6],"ko",'markersize',10);
xx=lsode("eom_w",[1.4-pi;6],tt);
plot(xx(:,1),xx(:,2),"k-",'linewidth',5)
xx=lsode("eom_w",[1.6-pi;6],tt);
plot(xx(:,1),xx(:,2),"k-",'linewidth',5)
xlabel('w=x-{\pi}','fontsize',42);
ylabel("dw/dt",'fontsize',42); drawnow();

Code 5 “expAt.m”

# "expAt.m"
A1=[0,1;-9.8,-1]; A2=[0,1;9.8,-1];
tt=linspace(0,10,200); xx=zeros(2,200);
x0=[0.5;-4.8];
for i=1:200
 xx(:,i)=expm(A1*tt(i))*x0;
endfor
clf; hold off;
plot(xx(1,:),xx(2,:),"k-",'linewidth',5);
axis([-2,8,-6,6]);
hold on;
x0a=[1.4-pi;6]; x0b=[1.6-pi;6];
for i=1:200
 xxa(:,i)=expm(A2*tt(i))*x0a;
 xxb(:,i)=expm(A2*tt(i))*x0b;
endfor
plot(xxa(1,:).+pi,xxa(2,:),"k-",'linewidth',5);
plot(xxb(1,:).+pi,xxb(2,:),"k-",'linewidth',5);
xlabel("x",'fontsize',32); ylabel("dx/dt",'fontsize',32); drawnow();

Code 6 “expAt2.m”

# "expAt2.m"
global A;
A=[0,1;-9.8,-1];
function dX=eqn(X,t)
 global A;
 Xmat=reshape(X,size(A)); #行列化
 dXmat=A*Xmat; #行列微分方程式
 dX=dXmat(:); #ベクトル化
endfunction
t=1; n=100; tt=linspace(0,t,n);
X0mat=eye(size(A));
X0=X0mat(:); #ベクトル化
XX=lsode("eqn",X0,tt);
Xmat=reshape(XX(n,:),size(A)) #手製の数値解
expm(A) #組み込みの行列指数関数

Code 7 “modexp.m”

# "modexp.m"
global A Adiag;
A=[0,1;-2,-3];
v1=[1;-1]; v1=v1/norm(v1); #unit vector
v2=[1;-2]; v2=v2/norm(v2); #unit vector
T=[v1,v2]; #change of basis matrix
function dx=eos(x,t)
 global A;
 dx = A*x;
endfunction
Adiag=inv(T)*A*T;
function dy=eos_mode(y,t)
 global Adiag;
 dy = Adiag*y;
endfunction
x0=[1;0];
y0=inv(T)*x0;
n=100;
tt=linspace(0,10,n);
xx=lsode("eos", x0, tt);
yy=lsode("eos_mode", y0, tt);
zz=zeros(n,2);
for i=1:n
 zz(i,:)=(yy(i,1)*v1 + yy(i,2)*v2)';
endfor
clf; subplot(1,2,1);
plot( xx(:,1), xx(:,2), "k-;x(t);" ); hold on;
plot( zz(:,1), zz(:,2), "ko;y_1(t)*v_1+y_2(t)*v_2;" );
xlabel("x_1(t)"); ylabel("x_2(t)");
subplot(1,2,2);
plot( yy(:,1), yy(:,2), "k-;y(t);" );
xlabel("y_1(t)"); ylabel("y_2(t)"); drawnow();

Code 8 “Wmat.m”

※Code 9 の実行に必要．

# "Wmat.m"
function W = Wmat( aa )
  n = length(aa);  #固有方程式の次数
  W=zeros(n,n);
  for i=1:n
    for j=1:n
      if ( i+j < n+1 )
        W(i,j) = aa(i+j);
      elseif ( i+j == n+1 )
        W(i,j) = 1;
      endif
    endfor
  endfor
endfunction

Code 9 “ctrbform.m”

※同じフォルダに，Code 8を書き込んだファイル Wmat.m が必要．

※実行結果の -0 は 0 のこと．

# "ctrbform.m"
function [y] = chop ( x ) #計算機誤差の除去
  eps=1e-12; y=( abs(x) > eps).*x;
endfunction
### システム行列 ###
n=3;
A=randn(n,n)
bb=randn(n,1)
### 可制御性行列 ###
Uc=[bb,A*bb,A*A*bb];
### 行列Aの固有方程式 ###
aa=poly(A); # aa(1)*s^n +aa(2)*s^(n-1) +... +aa(n)*s +aa(n+1), a(1)=1
aa=fliplr(aa)(1:n) #テキスト順 s^n +aa(n)*s^(n-1) +... +aa(1)
W = Wmat(aa)
### 基底変換行列 ###
Tc=Uc*W
### コンパニオン行列 ###
Ac=chop(inv(Tc)*A*Tc)
bbc=chop(inv(Tc)*bb)

Code 10 “eigplace.m”

# "eigplace.m"
global A bb KKp
### システム行列 ###
n=3;
A=randn(n,n);
bb=randn(n,1);
### 可制御性行列 ###
Uc=[bb,A*bb,A*A*bb];
### 行列Aの固有方程式 ###
aa=poly(A); # aa(1)*s^n +aa(2)*s^(n-1) +... +aa(n)*s +aa(n+1), a(1)=1
aa=fliplr(aa)(1:n); #テキスト順 s^n +aa(n)*s^(n-1) +... +aa(1)
W = Wmat(aa)
### 基底変換行列 ###
Tc=Uc*W;
### 目標の固有値と固有方程式 ###
ss=[-1-2i,-1+2i,-3]
cc=poly(diag(ss)); cc=fliplr(cc)(1:n); #テキスト順
### ゲイン ###
KKp=(cc-aa)*inv(Tc);
### 検算(閉ループ系の固有値) ###
eig(A-bb*KKp)
### シミュレーション ###
function dx=eos(x,t)
 global A bb KKp;
 ut = -KKp*x;
 dx = A*x + bb*ut;
endfunction
tt=linspace(0,10,200);
xx=lsode("eos",[1;1;1],tt);
plot(tt,xx,[";x1(t);";";x2(t);";";x3(t);"]);
xlabel("Time"); drawnow;

Code 11 “mck.m”

# "mck.m"
global k c omeg;
k=1; c=0.1;
### シミュレーションから計測した振幅比 ###
function dx=eom(x,t)
 global k c omeg;
 A=[ 0,  1; ...
    -k, -c ];
 b=[0;1];
 ut=cos(omeg*t);
 dx=A*x+b*ut;
endfunction
omn=51; freq1=linspace(0.6,1.4,omn);
xpp=[];
for omi=1:omn;
 omeg=freq1(omi);
 T=2*pi/omeg;
 tt=linspace(0,20*T,20*50);
 xx=lsode(@eom,[0;0],tt);
 tti=[20*45:20*50]; #終端時刻付近のデータ番号
 steady=xx(tti,1);  #定常応答
 plot(tt,cos(omeg*tt), "-", ... #入力
      tt,xx(:,1),"-",'linewidth',4,  ... #応答
      tt(tti),steady,"o"); #p-pの計測点
 xlabel("Time"); ylabel("x(t)");
 axis([0,20*T,-10,10]);
 title(sprintf("{/Symbol w}=%f",omeg));
 drawnow(); sleep(0.1);
 switch (omi)
  case 1
   title(""); print -deps2 -F:24 mck_01.eps
  case 26
   title(""); print -deps2 -F:24 mck_26.eps
  case 51
   title(""); print -deps2 -F:24 mck_51.eps
 endswitch
 pp=max(steady)-min(steady);
 xpp=[xpp;pp/2];
endfor
sleep(1);
plot(freq1,xpp,"o;x_{p-p}/2;",'markersize',8);
xlabel("{/Symbol w}"); ylabel("R"); axis([0.6,1.4,0,12]);
sleep(2);
### 伝達関数で計算した振幅比 ###
function y=G(s)
 global k c;
 A=[ 0,  1; ...
    -k, -c ];
 b=[0;1];
 I=eye(2); #単位行列
 inv(s*I-A);
 y=inv(s*I-A)*b;
endfunction
freq2=linspace(0.6,1.4,100); RR=[];
for j=1:100
  omeg=freq2(j);
  Giw=G(i*omeg);
  R=abs(Giw);   #振幅比 (変位,速度)
  phi=arg(Giw); #位相差 (変位,速度)
  RR=[RR,R(1)]; #振幅比 (変位)
endfor
hold on;
plot(freq2,RR,"r-;|G(i{/Symbol w})|;",'linewidth',5);
xlabel('{/Symbol w}'); ylabel('R'); axis([0.6,1.4,0,12]);
drawnow();
print -deps2 -F:24 mck_R.eps

Code 12 “dyndamp.m”

# "dyndamp.m"
global k c kmu cmu mu;
k=1; c=0.01; kmu=1; cmu=0.01;
mu=0.05;
function y=G(s)
 global k c kmu cmu mu;
 A=[      0,      1,      0,      0; ...
     -k-kmu, -c-cmu,     kmu,     cmu; ...
          0,      0,      0,      1; ...
      kmu/mu,  cmu/mu, -kmu/mu, -cmu/mu ];
 b=[0;1;0;0];
 I=eye(4); #単位行列
 inv(s*I-A);
 y=inv(s*I-A)*b;
endfunction
### 2つ目のばねが硬い場合 ###
kmu=100; #硬くして一体で動かす
omm=linspace(0.5,1.5,100);
RR1=[]; Phi1=[];
for j=1:100
  omeg=omm(j);
  Giw=G(i*omeg);
  Giw(3)=Giw(3)-Giw(1); #相対変位
  Giw(4)=Giw(4)-Giw(2); #相対速度
  R=abs(Giw)';   #振幅比
  phi=arg(Giw)'; #位相差
  RR1=[RR1;R]; Phi1=[Phi1;phi];
endfor
### 2つ目のばね・減衰器を最適調整した場合 ###
kmu=mu/(1+mu)^2
cmu=kmu*sqrt(3/2*mu*(1+mu))
RR2=[]; Phi2=[];
for j=1:100
  omeg=omm(j);
  Giw=G(i*omeg);
  phi=arg(Giw)'; #位相差
  Giw(3)=Giw(3)-Giw(1); #相対変位
  Giw(4)=Giw(4)-Giw(2); #相対速度
  R=abs(Giw)';   #振幅比 ↑ ここで差をとってはダメ
  RR2=[RR2;R]; Phi2=[Phi2;phi];
endfor
### 振幅比のグラフ ###
plot(omm,RR1(:,1),"r-;R_1;",'linewidth',5); hold on;
plot(omm,RR1(:,3),"bo;R'_3;",'markersize',5); hold off;
xlabel('{/Symbol w}'); ylabel('R'); axis([0.5,1.5,-4,120]);
drawnow(); sleep(1);
print -deps2 -F:26 dyndamp1.eps
plot(omm,RR2(:,1),"r-;R_1;",'linewidth',5); hold on;
plot(omm,RR2(:,3),"bo;R'_3;",'markersize',5); hold off;
xlabel('{/Symbol w}'); ylabel('R'); #axis([0.5,1.5,-1,25]);
drawnow(); sleep(1);
print -deps2 -F:26 dyndamp2.eps
### 位相差のグラフ ###
plot(omm,Phi1(:,1),"r-;{/Symbol f}_1;",'linewidth',5); hold on;
Phi2(:,3)=mod(Phi2(:,3).- 2*pi,2*pi)-2*pi;
plot(omm,Phi1(:,3),"bo;{/Symbol f}'_3;",'markersize',5); hold off;
xlabel('{/Symbol w}'); ylabel('{/Symbol f}'); axis([0.5,1.5,-pi,0]);
grid on; drawnow(); sleep(1);
print -deps2 -F:26 dyndamp_phi1.eps
plot(omm,Phi2(:,1),"r-;{/Symbol f}_1;",'linewidth',5); hold on;
Phi2(:,3)=mod(Phi2(:,3).- 2*pi,2*pi)-2*pi;
plot(omm,Phi2(:,3),"bo;{/Symbol f}'_3;",'markersize',5); hold off;
xlabel('{/Symbol w}'); ylabel('{/Symbol f}'); axis([0.5,1.5,-2*pi,0]);
grid on; drawnow();
print -deps2 -F:26 dyndamp_phi2.eps
save dyndamp.dat RR1 RR2 Phi2

Code 13 “optbvp.m”

# "optbvp.m"
a=1; b=1; t1=1; q=100; r=1;
global A;
function dx=eom(x,t)
 global A;
 dx = A*x;
endfunction
A=[a, -b^2/(2*r); -2*q, -a];
expA=expm(A*t1);
aa1=expA(:,1); aa2=expA(:,2); ee1=[1;0]; ee2=[0;1];
tt=linspace(0,t1,200);
# 両端固定 ... x0 x1 既知
x0=1; x1=2; p0p1=[-aa2,ee2]\(x0*aa1-x1*ee1); ## p ... lamda
xx=lsode("eom",[x0;p0p1(1)],tt);
uu=-b/(2*r)*xx(:,2);
subplot(2,2,1); plot(tt,xx(:,1),"k-;x(t);");
title("Both sides fixed",'fontsize',24);
subplot(2,2,3); plot(tt,uu,"k-;u(t);");
xlabel("Time",'fontsize',24);
# 終端自由 ... x0 p1 既知
x0=1; p1=0; p0x1=[-aa2,ee1]\(x0*aa1-p1*ee2);
xx=lsode("eom",[x0;p0x1(1)],tt);
uu=-b/(2*r)*xx(:,2);
subplot(2,2,2); plot(tt,xx(:,1),"k-;x(t);");
title("Free endpoint",'fontsize',24);
subplot(2,2,4); plot(tt,uu,"k-;u(t);");
xlabel("Time",'fontsize',24);
xxb=xx(:,1); uub=uu; ttb=tt; save "optbvp.dat" xxb uub ttb;

Code 14 “optRic.m”

※先に，同じフォルダで Code 13 を実行しておく必要がある．

# "optRic.m"
global a b q r tt KK;
a=1; b=1; q=100; r=1; x0=1; t1=1; #t1=5;
function dx=eos(x,t) #状態方程式
 global a b;
 u = - K(t)*x;
 dx = a*x + b*u;
endfunction
function dP=Riccati(P,t) #リッカチ方程式
 global a b q r;
 dP = -q - a*P - P*a + P*b/r*b*P;
endfunction
function y = K(t) #KK(1),...,KK(n) を補間して K(t) を求める
 global tt KK;
 y = interp1(tt, KK, t);
endfunction
n=100; tt=linspace(0,t1,n);
# リッカチ方程式を解く
rt=tt(n:-1:1);               #逆時間
PP=lsode("Riccati",0,rt);    #P(t) の数列
KK=1/r*b*PP(n:-1:1);         #K(t) の数列
clf; subplot(3,1,1); plot(rt, PP(:,1), "k-;;", 'linewidth', 1);
ylabel("P(t)",'fontsize',24);
# 状態方程式を解く (終端自由)
xx=lsode("eos",x0,tt);       #x(t) の数列
load "optbvp.dat";           #optbvp.m のデータを読む
subplot(3,1,2);
plot(tt,xx,"o;with P(t);",'markersize',8,...
     ttb,xxb,"k-;with lam(t);", 'linewidth', 1);
ylabel("x(t)",'fontsize',24);
subplot(3,1,3); uu=-KK.*xx;
plot(tt,uu,"o;with P(t);",'markersize',8,...
     ttb,uub,"k-;with lam(t);", 'linewidth', 1);
xlabel("t",'fontsize',24); ylabel("u(t)",'fontsize',24);

Code 15 “optlqr.m”

# "optlqr.m"
global A bb KK
### システム行列 ###
n=3;
A=randn(n,n);
bb=randn(n,1);
### 無限次元最適レギュレータ ###
Q=diag([1,1,1]); #状態x(t)の重み
R=2;             #入力u(t)の重み
[KK,PP,EE]=lqr(A,bb,Q,R) #最適レギュレータを求める
### シミュレーション ###
function dx=eos(x,t)
 global A bb KK;
 ut = -KK*x;
 dx = A*x + bb*ut;
endfunction
tt=linspace(0,10,200);
xx=lsode("eos",[1;1;1],tt);
plot(tt,xx,["-0;x1(t);";"-1;x2(t);";"-3;x3(t);"],'linewidth',2);
xlabel("Time"); drawnow;

Code 16 “invp.m”

# "invp.m"
global M m l g KK;  #大域変数の宣言
M=2/3; m=1/3; l=1; g=9.8;
KK=[0,0,0,0]; #制御無し
function dx=eqn(x,t)
 global M m l g KK; #関数の中から大域変数を参照する宣言
 ut = -KK*x;  #制御入力
 D=M*l + m*l*(1-cos(x(3))^2);
 dx(1) = x(2);
 dx(2) = ( m*l^2*x(4)^2*sin(x(3)) - m*l*g*cos(x(3))*sin(x(3)) )/D ...
         + l/D * ut;
 dx(3) = x(4);
 dx(4) = ( -m*l^2*x(4)^2*cos(x(3))*sin(x(3)) + (M+m)*g*sin(x(3)) )/D ...
         - cos(x(3))/D * ut;
endfunction
x0=[0; 0; 0.2; 0];       #初期値(ちょっと倒れて静止)
n=70;                    #時間のきざみ数
tt=linspace(0,7,n);     #時間0～10の等比数列
xx=lsode("eqn", x0, tt); #常微分方程式"eqn"を解く
for i=1:n  #簡易アニメーション
 x=xx(i,1);              #第1変数が台車の位置
 theta=xx(i,3);          #第3変数が倒れ角
 plot( [x,x+l*sin(theta)], [0,l*cos(theta)] ); #線分のプロット
 axis([-5,5,-5,5],'square');    #座標軸の設定
 title(sprintf("xxxxx: %d",i)); #タイトルの設定
 grid on; drawnow();
 if ( i==1 ) sleep(1);   #初回は1秒待機
 else sleep(0.1); endif  #この数値でアニメーション速度を調整
endfor

Code 17 “invpGain.m”

# "invpGain.m"
global A bb KKp
### システム行列 (倒立振子)###
n=4;
M=2/3; m=1/3; l=1; g=9.8;
A=[0,1,0,0; 0,0,-m/M*g,0; 0,0,0,1; 0,0,(M+m)/(M*l)*g,0];
bb=[0;1/M;0;-1/(M*l)];
### 可制御性行列 ###
Uc=[bb,A*bb,A*A*bb,A*A*A*bb]; #*** 4次元化 ***
### 行列Aの固有方程式 ###
aa=poly(A); # aa(1)*s^n +aa(2)*s^(n-1) +... +aa(n)*s +aa(n+1), a(1)=1
aa=fliplr(aa)(1:n); #テキスト順 s^n +aa(n)*s^(n-1) +... +aa(1)
W = Wmat(aa)
### 基底変換行列 ###
Tc=Uc*W;
### 目標の固有値と固有方程式 ###
ss=[-1,-2,-1-2*i,-1+2*i]  #*** 4次元化 ***
cc=poly(diag(ss)); cc=fliplr(cc)(1:n); #テキスト順
### ゲイン ###
KKp=(cc-aa)*inv(Tc)
### 検算(閉ループ系の固有値) ###
eig(A-bb*KKp)
### シミュレーション ###
function dx=eos(x,t)
 global A bb KKp;
 ut = -KKp*x;
 dx = A*x + bb*ut;
endfunction
tt=linspace(0,10,200);
xx=lsode("eos",[0;0;0.2;0],tt); #*** 4次元化 ***
plot(tt,xx,[";x1(t);";";x2(t);";";x3(t);";";x4(t);"]); #*** 4次元化 ***
xlabel("Time");drawnow;
save invpGain.dat KKp;

Code 18 “invpLQR.m”

# "invpLQR.m"
global A bb KKopt
### システム行列 (倒立振子)###
n=4;
M=2/3; m=1/3; l=1; g=9.8;
A=[0,1,0,0; 0,0,-m/M*g,0; 0,0,0,1; 0,0,(M+m)/(M*l)*g,0];
bb=[0;1/M;0;-1/(M*l)];
### 無限時間最適レギュレータ ###
Q=diag([10,3,7,3]); #状態x(t)の重み
R=2;                #入力u(t)の重み
[KKopt,PP,EE]=lqr(A,bb,Q,R) #最適レギュレータを求める
### シミュレーション ###
function dx=eos(x,t)
 global A bb KKopt;
 ut = -KKopt*x;
 dx = A*x + bb*ut;
endfunction
tt=linspace(0,10,200);
xx=lsode("eos",[0;0;0.2;0],tt);
plot(tt,xx,[";x1(t);";";x2(t);";";x3(t);";";x4(t);"]);
xlabel("Time");drawnow;
save invpLQR.dat KKopt;

Code 19 “pendvec-pf.m”

※同じフォルダに，Code 2を書き込んだファイル vectorfield.m が必要．

# "pendvec-pf.m"
global c k K;
c=1; k=9.8; K=20;
function dx=eom(x,t)
 global c k K;
 u = -K*x(1);
 dx(1) = x(2);
 dx(2) = -c*x(2) + k*sin(x(1)) + u;
endfunction
function myplot(nx, ny, rx, ry, x0s)
  vectorfield("eom",nx,ny,rx,ry,1.2); hold on;
  tt=linspace(0,30,500);
  XX=[]; YY=[];
  for i=1:columns(x0s)
   xx=lsode("eom",x0s(:,i),tt);
   XX=[XX,xx(:,1)]; YY=[YY,xx(:,2)];
  endfor
  plot(XX,YY,"k-",'linewidth',2)
  axis([rx(1),rx(2),ry(1),ry(2)]);
  xlabel("x"); ylabel("dx/dt");
endfunction
clf; hold off;
subplot(2,2,1);
K=20; rx=[-2,2]; ry=[-6,6];
x01=[-1.5;-4]; x02=[1.6;3];
myplot(20,15,rx,ry,[x01,x02]); drawnow();
subplot(2,2,2);
K=9.75; rx=[-2,2]; ry=[-6,6];
x01=[-1.5;-2]; x02=[1.5;2];
myplot(20,15,rx,ry,[x01,x02]);
plot(0.15*[0;1;-1],[0;0;0],"o",'markersize',5);
drawnow();
subplot(2,2,3);
K=8; ry=[-3,3];
x01=[-1.8;1.5]; x02=[-1.8;2];
myplot(20,15,rx,ry,[x01,x02]);
plot(1.08*[0;1;-1],[0;0;0],"o",'markersize',5);

Code 20 “pendvec-sn.m”

※同じフォルダに，Code 2を書き込んだファイル vectorfield.m が必要．

# "pendvec-sn.m"
global c k K F;
c=1; k=9.8; K=8;
function dx=eom(x,t)
 global c k K F;
 u = -K*x(1) + F*cos(x(1));
 dx(1) = x(2);
 dx(2) = -c*x(2) + k*sin(x(1)) + u;
endfunction
function myplot(nx, ny, rx, ry, x0s)
  vectorfield("eom",nx,ny,rx,ry,1.2); hold on;
  tt=linspace(0,30,500);
  XX=[]; YY=[];
  for i=1:columns(x0s)
   xx=lsode("eom",x0s(:,i),tt);
   XX=[XX,xx(:,1)]; YY=[YY,xx(:,2)];
  endfor
  plot(XX,YY,"k-",'linewidth',2)
  axis([rx(1),rx(2),ry(1),ry(2)]);
  xlabel("x"); ylabel("dx/dt");
endfunction
clf; hold off;
subplot(2,2,1);
F=0.9; rx=[-1.5,1.5];ry=[-3,3];
x01=[0.53;-2]; x02=[0.5;-2];
myplot(20,15,rx,ry,[x01,x02]); drawnow();
plot([-0.58;-0.81],[0;0],"o",'markersize',8); drawnow();
subplot(2,2,2);
F=2; rx=[-2,2];ry=[-3,3];
x01=[-1.4;-1]; x02=[0.235;-2];
myplot(20,15,rx,ry,[x01,x02]); drawnow();
subplot(2,2,3);
F=-0.9; rx=[-1.5,1.5];ry=[-3,3];
x01=-[0.53;-2]; x02=-[0.5;-2];
myplot(20,15,rx,ry,[x01,x02]); drawnow();
plot(-[-0.58;-0.81],[0;0],"o",'markersize',8); drawnow();
subplot(2,2,4);
F=-2; rx=[-2,2];ry=[-3,3];
x01=-[-1.4;-1]; x02=-[0.235;-2];
myplot(20,15,rx,ry,[x01,x02]); drawnow();

Code 21 “pend-jump.m”

# "pend-jump.m"
global c k K FR dF;
c=0.7; k=9.8;
function dx=eom(x,t)
 global c k K FR dF;
 u = -K*x(1);
 F = FR(1)+dF*t;
 dx(1) = x(2);
 dx(2) = -c*x(2) + k*sin(x(1)) + F*cos(x(1)) + u;
endfunction
K=9.5; T=200; n=120;
x0=[-0.4;0]; FR=[-0.3,0.3]; dF=(FR(2)-FR(1))/T;
tt=linspace(0,T,n);
x1=lsode("eom",x0,tt);
F1=FR(1) .+ tt*dF;
x0=x1(n,:)'; FR=-FR; dF=(FR(2)-FR(1))/T;
x2=lsode("eom",x0,tt);
F2=FR(1) .+ tt*dF;
xx=[x1;x2]; FF=[F1,F2];
for i=1:2*n
  x=sin(xx(i,1)); y=cos(xx(i,1));
  plot([0,x],[0,y],"-",'linewidth',2);
  axis([-1.5,1.5,-1.5,1.5],'square');
  title(sprintf("F=%f",FF(i)));  drawnow();
  if ( i==1 ) sleep(2);
  else sleep(0.02); endif
endfor
plot(F1,x1(:,1),"-;F increases -->;",'linewidth',2, ...
     F2,x2(:,1),"--;F decreases <--;",'linewidth',2);
xlabel("F"); ylabel("x");
axis([-0.3,0.3,-1,1]); drawnow();

Code 22 “crane-hopf.m”

# "crane-hopf.m" ※旧"invp-hopf.m"
global M m l g KK;  #大域変数の宣言
M=2/3; m=1/3; l=1; g=-9.8;
KK=[22.4933,7,3,5]; #Hopf分岐点
KK(1)=0.8*KK(1) #リミットサイクル側にちょっとずらす
function dx=eqn(x,t)
 global M m l g KK; #関数の中から大域変数を参照する宣言
 ut = -KK*x;  #制御入力
 D=M*l + m*l*(1-cos(x(3))^2);
 dx(1) = x(2);
 dx(2) = ( m*l^2*x(4)^2*sin(x(3)) - m*l*g*cos(x(3))*sin(x(3)) )/D ...
         + l/D * ut;
 dx(3) = x(4);
 dx(4) = ( -m*l^2*x(4)^2*cos(x(3))*sin(x(3)) + (M+m)*g*sin(x(3)) )/D ...
         - cos(x(3))/D * ut;
endfunction
x0=[1; 0.2; 0; 0];       #初期値(ちょっと倒れて静止)
n=1500;                   #時間のきざみ数
tt=linspace(0,30,n);     #時間0〜10の等比数列
xx=lsode("eqn", x0, tt); #常微分方程式"eqn"を解く
subplot(2,1,1); plot( xx(:,1), xx(:,2) );
subplot(2,1,2); plot( xx(:,3), xx(:,4) );
sleep(2);clf;
for i=1:5:n  #簡易アニメーション 5つ飛ばしで描画
 x=xx(i,1);              #第1変数が台車の位置
 theta=xx(i,3);          #第3変数が倒れ角
 plot( [x,x+l*sin(theta)], [0,l*cos(theta)] ); #線分のプロット
 axis([-5,5,-5,5],'square');    #座標軸の設定
 title(sprintf("xxxxx: %d",i)); #タイトルの設定
 grid on; drawnow();
 if ( i==1 ) sleep(1);   #初回は1秒待機
 else sleep(0.02); endif  #この数値でアニメーション速度を調整
endfor

Code 23 “newton.m”

# "newton.m"
function y=f(x)
 y=x^2-1;
endfunction
function y=df(x)
 y=2*x;
endfunction
function xx=newton(x0,n)
 x=x0;
 xx=[x];
 for i=1:n
  x = x - f(x)/df(x);
  xx = [xx,x];
 endfor
endfunction
hold on; rand('seed',1);
for i=1:50
 x0=10*rand(1)-5; # -5<x<5 の一様乱数
 xx=newton(x0,10);
 plot(0:10,xx);
 xlabel('n'); ylabel('x_n');
endfor
#print -deps -F:24 newton.eps;

Code 24 “circle.m”

# circle.m
function z=f(x,y)
 z=x^2+y^2-1;
endfunction
dx=0.05; dy=0.05;
x=-1/sqrt(2); y=-x;
xx=[x];yy=[y];
while (x<0.8)
 x=x+dx; y=fsolve(@(y) f(x,y),y); #1変数関数 f(y) にする文法
 xx=[xx,x]; yy=[yy,y];
endwhile
hold on;
plot(xx,yy,"o");
xx=[x];yy=[y];
while (y>-0.8)
 y=y-dy; x=fsolve(@(x) f(x,y),x); #1変数関数 f(x) にする文法
 xx=[xx,x]; yy=[yy,y];
endwhile
plot(xx,yy,"x");
xx=[x];yy=[y];
while (x>-0.8)
 x=x-dx; y=fsolve(@(y) f(x,y),y);
 xx=[xx,x]; yy=[yy,y];
endwhile
plot(xx,yy,"*");
xx=[x];yy=[y];
while (y<0.8)
 y=y+dy; x=fsolve(@(x) f(x,y),x);
 xx=[xx,x]; yy=[yy,y];
endwhile
plot(xx,yy,"+"); xlabel("x"); ylabel("y");
#print -deps -F:24 circle.eps

Code 25 “bdsn-ini.m”

# "bdsn-ini.m"
global c k K F;
c=0.7; k=9.8; K=20;
function dx=eom(x,t)
 global c k K F;
 u = -K*x(1);
 dx(1) = x(2);
 dx(2) = -c*x(2) + k*sin(x(1)) + F*cos(x(1)) + u;
endfunction
format long; K=9.5; F0=-0.3
### 平衡点の概算 ###
tt=linspace(0,200,120);
F=F0; xx=lsode("eom",[0;0],tt);
subplot(2,1,1); plot(tt,xx(:,1)); xlabel("t"); ylabel("x1");
subplot(2,1,2); plot(tt,xx(:,2)); xlabel("t"); ylabel("x2");
drawnow(); sleep(2);
x0a=xx(120,1) #最後の値=概算値
### 初期解 ###
function y=f(x)
 global c k K F;
 u = -K*x;
 y = k*sin(x) + F*cos(x) + u;
endfunction
x0=fsolve("f",x0a) #高精度化したもの
save bdsn-ini.dat c k K F0 x0 x0a;

Code 26 “bdsn-cont.m”

# "bdsn-cont.m"
global c k K;
load bdsn-ini.dat; #c k K F0 x0 x00;
function y=f(x,F)
 global c k K;
 u = -K*x;
 y = k*sin(x) + F*cos(x) + u;
endfunction
### 解の接続 ###
dx=0.03; dF=0.02;
x=x0; F=F0;
xx=[x]; FF=[F];
while (F<0.04) #0.04はCode 19のグラフのカーブより目測で
 F=F+dF; x=fsolve(@(x) f(x,F),x); #掃引方向→
                  #1変数関数 f(x) にする文法
 xx=[xx,x]; FF=[FF,F];
endwhile
hold on; plot(F,x,"ro",'markersize',8); #掃引方向の切替点
while (x<0.4) #0.4はこのプログラムの実行結果を見ながら調整
 x=x+dx; F=fsolve(@(F) f(x,F),F); #掃引方向↑
                  #1変数関数 f(F) にする文法
 xx=[xx,x]; FF=[FF,F];
endwhile
plot(F,x,"ro",'markersize',8); #掃引方向の切替点
while (F<0.3) #0.3はCode 19のグラフに合わせて
 F=F+dF; x=fsolve(@(x) f(x,F),x); #掃引方向→
 xx=[xx,x]; FF=[FF,F];
endwhile
### 分岐図 ###
plot(FF,xx,"-",'linewidth',3); #分岐図
xlabel("F"); ylabel("x");
axis([-0.3,0.3,-1,1]); drawnow;
save bdsn-cont.dat c k K FF xx
#print -mono -deps -F:24 bdsn-cont.eps

Code 27 “bdsn-stab.m”

# "bdsn-stab.m"
global c k K;
load bdsn-cont.dat #c k K FF xx
function m=Jaco(x,F) #ヤコビ行列
 global c k K;
 m = [0, 1; ...
      k*cos(x)-F*sin(x)-K, -c];
endfunction
function y=eigmax(x,F) #最大の固有値実部
 y=max(real(eig(Jaco(x,F))));
endfunction
### 安定判別 ###
xs=[]; Fs=[]; xu=[]; Fu=[];
for i=1:length(xx);
 stab=eigmax(xx(i),FF(i));
 if (stab<0)
  xs=[xs,xx(i)]; Fs=[Fs,FF(i)]; #安定平衡点
 else
  xu=[xu,xx(i)]; Fu=[Fu,FF(i)]; #不安定平衡点
 endif
endfor
### プロット用に安定平衡点を2分割する ###
jump=0; ns=length(xs);
for i=1:ns-1;
 if (abs(xs(i+1)-xs(i))>0.5)
  jump=i;
 endif
endfor
xs1=xs(1:jump); Fs1=Fs(1:jump);
xs2=xs(jump+1:ns); Fs2=Fs(jump+1:ns);
### 分岐図の描画 ###
hold on;
plot(Fs1,xs1,"-",'linewidth',3); #安定
plot(Fs2,xs2,"-",'linewidth',3); #安定
plot(Fu,xu,".1",'linewidth',3);  #不安定
xlabel("F"); ylabel("x");
axis([-0.3,0.3,-1,1]); drawnow;
save bdsn-stab.dat c k K Fs1 xs1 Fs2 xs2;
#print -deps2 -F:26 bdsn-stab.eps

Code 28 “crane-hopfbd.m”

# "crane-hopfbd.m" ※旧"invp-hopfbd.m"
global A bb K3 K4;
K3=3; K4=5; #台車制御のゲイン(固定)
M=2/3; m=1/3; l=1; g=-9.8;
A=[0,1,0,0; 0,0,-m/M*g,0; 0,0,0,1; 0,0,(M+m)/(M*l)*g,0];
bb=[0;1/M;0;-1/(M*l)];
### 閉ループ系の行列 ###
function m=H(K1,K2)
 global A bb K3 K4;
 KK=[K1,K2,K3,K4];
 m=A-bb*KK; #
endfunction
### 分岐方程式 ###
function y=bifeqn(K1,K2)
 m=H(K1,K2); #閉ループ系
 aa=poly(m)(2:rows(m)+1);
 y=aa(3)^2-aa(1)*aa(2)*aa(3)+aa(1)^2*aa(4);
endfunction
### Hopf分岐点の概算 ###
K2=7; KK1=linspace(20,25,100); #invp-hopf.mを参考に
yy=[];
for i=1:100
 yy=[yy,bifeqn(KK1(i),K2)];
endfor
plot(KK1,yy); xlabel("K1"); ylabel("bifeqn");
grid on; sleep(3); #グラフを見ると22.5付近でy=0
### 初期解 ###
format long;
K1=fsolve(@(K1) bifeqn(K1,K2),22.5) #概算より
          #1変数関数 bifeqn(K1) にする文法
BF1=[K1]; BF2=[K2];
### 接続1 ###
dK=0.1;
while (K2<=30)
 K2=K2+dK; K1=fsolve(@(K1) bifeqn(K1,K2),K1);
 BF1=[BF1,K1]; BF2=[BF2,K2];
endwhile
### 接続2 ###
K1=BF1(1); K2=BF2(1);
BF1=fliplr(BF1); BF2=fliplr(BF2); #接続1とグラフを継げたいので，逆順に
while (K1<=30)
 K1=K1+dK; K2=fsolve(@(K2) bifeqn(K1,K2),K2);
 BF1=[BF1,K1]; BF2=[BF2,K2];
endwhile
plot(BF1,BF2,"k-",'linewidth',5);
axis([5,30,5,30]);
xlabel("K1"); ylabel("K2");
### 固有値の確認 ###
hold on; format short;
p0=[BF1(170),BF2(170)]; #分岐点の1つ
ps=p0+[2,2]; pu=p0-[2,2];
plot(ps(1),ps(2),"k*",'markersize',10);
plot(p0(1),p0(2),"ko",'markersize',10);
plot(pu(1),pu(2),"k+",'markersize',10);
printf("(*)"); Stable_eigenvalue=eig(H(ps(1),ps(2)))
printf("(o)"); Imaginary_eigenvalue=eig(H(p0(1),p0(2)))
printf("(+)"); Unstable_eigenvalue=eig(H(pu(1),pu(2)))
#print -deps -F:24 invp-hopfbd.eps

---

以下，Octave の開発環境に含まれるコマンド mkoctfile を使う．

参考： http://www.kkaneko.com/rinkou/octave/octavemex.html

※例えば，Ubuntuでは「octave3.0-headers」みたいなパッケージを追加インストールすると使用できる．

Code 29 “fpen_po.cc”

Octave の開発環境で，Linux の場合，

> mkoctfile fpen_po.cc

として fpen_po.oct を作っておく．後述のプログラムで使う．※ 単独では使えない．

/* "fpen_po.cc"
 *
 * [...]$ mkoctfile fpen_po.cc
 *
 * として fpen_po.oct に変換する．
 */
#include <octave/oct.h>
#include <octave/LSODE.h>

/* 運動方程式
 * ※ octave でベクトルの第1成分は x(1) だが，ここでは x(0)
 * --- ここから --- */
struct _p { double c, k, P, omeg; } p;
ColumnVector xdot(const ColumnVector& x, double t) {
  ColumnVector dx(2);
  dx(0) = x(1); /* octave では dx(1)=x(2) */
  dx(1) = -p.c*x(1) -p.k*sin(x(0)) + p.P*cos(p.omeg*t);
  return dx;
}
/* --- ここまで --- */

DEFUN_DLD (fpen_po, args, , "Poincare map for a single pendlumn") {
  /* 引数の取得 */
  double get_global_double( char* s);
  ColumnVector x0 (args(0).vector_value ()); //初期値
  double T = (args(1).double_value ()); //周期
  double m = (args(2).int_value ()); //無視する過渡の点数
  double n = (args(3).int_value ()); //出力点数
  /* global 変数の取得 */
  p.c = get_global_double( "c" );
  p.k = get_global_double( "k" );
  p.P = get_global_double( "P" );
  p.omeg = get_global_double( "omeg" );
  /* 数値積分 */
  int dim = x0.rows(); //次元
  int Tn = 200; //分割数
  double t=0.0, dt = T/(double)Tn;
  ODEFunc odef(xdot);
  ColumnVector x=x0; //解
  for ( int i=1; i<m; i++ ) { /* 過渡を無視する */
    LSODE ls(x, 0.0, odef);
    for ( int ti=1; ti<=Tn; ti++ ) {
      t = dt*(double)ti;
      x = ls.do_integrate(t);
    }
  }
  Matrix pm(n,dim);
  for ( int j=0; j<dim; j++ ) {
    pm(0,j)=x(j);
  }
  for ( int i=1; i<n; i++ ) {
    LSODE ls(x, 0.0, odef);
    for ( int ti=1; ti<=Tn; ti++ ) {
      t = dt*(double)ti;
      x = ls.do_integrate(t);
    }
    for ( int j=0; j<dim; j++ ) {
      pm(i,j)=x(j);
    }
  }
  return octave_value (pm);
}

double get_global_double( char* s ) {
  double retval = 0.123456789012345;
  octave_value tmp = get_global_value (s, true);
  if ( ! tmp.is_defined () || tmp.is_empty() )
    octave_stdout<< "Global variable \"" << s << "\" is empty";
  else
    retval = tmp.double_value ();
  return retval;
}

Code 30 “fpen-traj.m”

# "fpend-traj.m" ※環境によっては実行時間がかかる！
global c k P omeg;
function dx=eom(x,t)
 global c k P omeg;
 dx(1) = x(2);
 dx(2) = -c*x(2) -k*sin(x(1)) + P*cos(omeg*t);
endfunction
k=1; P=2.7; omeg=1; T=2*pi/omeg;
x0=[0;0];
tn=1000; #無視する過渡応答
dn=tn+500; #取得するデータ長
tt=linspace(0,dn*T,50*dn);
### 1周期 ###
c=0.7; xx=lsode("eom",x0,tt)(tn+1:dn,:); #過渡を削除
hold off; plot(xx(:,1),xx(:,2),"-",'linewidth',5);
hold on; plot(xx(1,1),xx(1,2),"o",'markersize',10);
xlabel("x(t)"); ylabel('{~x{.5.}(t)}'); axis([-9,-2,-3,3]); drawnow();
print -deps2 -F:24 fpen-traj1.eps
### 2周期 ###
x0=xx(1,:)'; #分岐図と履歴を合わせるため
c=0.65; xx=lsode("eom",x0,tt)(tn+1:dn,:); #過渡を削除
hold off; plot(xx(:,1),xx(:,2),"-",'linewidth',5);
hold on; plot(xx([1,51],1),xx([1,51],2),"o",'markersize',10);
xlabel("x(t)"); ylabel('{~x{.5.}(t)}'); axis([-9,-2,-3,3]); drawnow();
print -deps2 -F:24 fpen-traj2.eps
### 4周期 ###
x0=xx(1,:)'; #分岐図と履歴を合わせるため
c=0.63; xx=lsode("eom",x0,tt)(tn+1:dn,:); #過渡を削除
hold off; plot(xx(:,1),xx(:,2),"-",'linewidth',5);
hold on; plot(xx([1,51,101,151],1),xx([1,51,101,151],2),"o",'markersize',10);
xlabel("x(t)"); ylabel('{~x{.5.}(t)}'); axis([-9,-2,-3,3]); drawnow();
print -deps2 -F:24 fpen-traj4.eps
### 無限周期 ###
x0=xx(1,:)'; #分岐図と履歴を合わせるため
c=0.61; xx=lsode("eom",x0,tt)(tn+1:dn,:); #過渡を削除
hold off; plot(xx(:,1),xx(:,2),"-",'linewidth',5);
xlabel("x(t)"); ylabel('{~x{.5.}(t)}'); axis([-9,-2,-3,3]); drawnow();
#print -deps2 -F:24 fpen-trajx.eps

Code 31 “fpen-bd.m”

※同じフォルダに，上で作った fpen_po.oct が必要．

# "fpen-bd.m"
global c k P omeg;
c=0.22; k=1; P=2.7; omeg=1.0;
n1=1000; #無視する過渡応答
n2=500;  #データ数
x0=[0;0]; N=100;
p=linspace(0.7,0.6,N+1);
hold on; fpen_bd=[];
for i=1:N+1
 c=p(i); T=2*pi/omeg;
 xx=fpen_po(x0,T,n1,n2); #fpen_po.oct によるポアンカレ写像の計算
 graph=[p(i)*ones(n2,1),xx(:,1)];
 plot(graph(:,1),graph(:,2),"."); drawnow();
 fpen_bd=[fpen_bd;graph];
 x0=xx(n2,:)';
endfor
xlabel('c', 'fontsize', 32); ylabel('x_n', 'fontsize', 32); drawnow();
save fpen-bd.dat fpen_bd;
#print -deps2 -F:26 fpen-bd.eps

Code 32 “fpen-poincare.m”

※同じフォルダに，上で作った fpen_po.oct が必要．

# "fpen-poincare.m"
global c k P omeg;
n1=1000; #無視する過渡応答
n2=4000; #データ数
x0=[0;0];
c=0.61; k=1; P=2.7; omeg=1.0;
T=2*pi/omeg; #外力の周期
xx=fpen_po(x0,T,n1,n2); #ポアンカレ写像
x=mod(xx(:,1).+pi,2*pi).-pi; #角度を-πからπに収める変換
y=xx(:,2);
plot(x,y,"."); axis([-pi,pi,-0.71,3.22]);
xlabel('x_n mod {2\pi}'); ylabel('{~x{.5.}_n}'); drawnow();
print -deps2 -F:26 fpen-poincare1.eps
c=0.3;
xx=fpen_po(x0,T,n1,n2); #ポアンカレ写像
x=mod(xx(:,1).+pi,2*pi).-pi; #角度を-πからπに収める変換
y=xx(:,2);
plot(x,y,"."); axis([-pi,pi,-0.71,3.22]);
xlabel('x_n mod {2\pi}'); ylabel('{~x{.5.}_n}'); drawnow();
#print -deps2 -F:26 fpen-poincare2.eps

Code 33 “fpen_polin.cc”

Octave の開発環境で，Linux の場合，

> mkoctfile fpen_polin.cc

として fpen_polin.oct を作っておく．後述のプログラムで使う．※ 単独では使えない．

/* "fpen_polin.cc"
 *
 * [...]$ mkoctfile fpen_polin.cc
 *
 * として fpen_polin.oct に変換する．
 */
#include <octave/oct.h>
#include <octave/LSODE.h>

/* 運動方程式
 * ※ octave でベクトルの第1成分は x(1) だが，ここでは x(0)
 * --- ここから --- */
struct _p { double c, k, P, omeg; } p;
ColumnVector xdot(const ColumnVector& x, double t) {
  ColumnVector dx(4);
  dx(0) = x(1); /* octave では dx(1)=x(2) */
  dx(1) = -p.c*x(1) -p.k*sin(x(0)) + p.P*cos(p.omeg*t);
  dx(2) = x(3); /* octave では dx(3)=x(4) */
  dx(3) = -p.k*cos(x(0))*x(2) -p.c*x(3) ;
  return dx;
}
/* --- ここまで --- */
/* ※ 以下は"fpen_po.cc"と共通 */

DEFUN_DLD (fpen_polin, args, , "Poincare map for a single pendlumn") {
  /* 引数の取得 */
  double get_global_double( char* s);
  ColumnVector x0 (args(0).vector_value ()); //初期値
  double T = (args(1).double_value ()); //周期
  double m = (args(2).int_value ()); //無視する過渡の点数
  double n = (args(3).int_value ()); //出力点数
  /* global 変数の取得 */
  p.c = get_global_double( "c" );
  p.k = get_global_double( "k" );
  p.P = get_global_double( "P" );
  p.omeg = get_global_double( "omeg" );
  /* 数値積分 */
  int dim = x0.rows(); //次元
  int Tn = 200; //分割数
  double t=0.0, dt = T/(double)Tn;
  ODEFunc odef(xdot);
  ColumnVector x=x0; //解
  for ( int i=1; i<m; i++ ) { /* 過渡を無視する */
    LSODE ls(x, 0.0, odef);
    for ( int ti=1; ti<=Tn; ti++ ) {
      t = dt*(double)ti;
      x = ls.do_integrate(t);
    }
  }
  Matrix pm(n,dim);
  for ( int j=0; j<dim; j++ ) {
    pm(0,j)=x(j);
  }
  for ( int i=1; i<n; i++ ) {
    LSODE ls(x, 0.0, odef);
    for ( int ti=1; ti<=Tn; ti++ ) {
      t = dt*(double)ti;
      x = ls.do_integrate(t);
    }
    for ( int j=0; j<dim; j++ ) {
      pm(i,j)=x(j);
    }
  }
  return octave_value (pm);
}

double get_global_double( char* s ) {
  double retval = 0.123456789012345;
  octave_value tmp = get_global_value (s, true);
  if ( ! tmp.is_defined () || tmp.is_empty() )
    octave_stdout<< "Global variable \"" << s << "\" is empty";
  else
    retval = tmp.double_value ();
  return retval;
}

Code 34 “fpen-stab.m”

※同じフォルダに，上で作った fpen_polin.oct が必要．

# "fpen-stab.m"
global c k P omeg T;
c=0.7; k=1; P=2.7; omeg=1.0;
n1=1000; #無視する過渡応答
n2=4000; #データ数
### 不動点の概算 ###
x0=[0;0];
T=2*pi/omeg; #外力の周期
xx=fpen_po(x0,T,n1,n2); #ポアンカレ写像
p0=xx(n2,:)'; #概算値
### 不動点の方程式 ###
function y=fixp(p0)
 global T;
 p1=fpen_po(p0,T,0,2)(2,:)';
 y=p0-p1;
endfunction
p=p0; nc=200; cc=linspace(0.7,0.6,nc+1);
fpen_stab=[];
for i=1:nc+1
 c=cc(i);
 p=fsolve("fixp",p); ### ニュートン法による不動点 ###
 R0=[p;1;0]; R1=fpen_polin(R0,T,0,2)(2,3:4)';
 R0=[p;0;1]; R2=fpen_polin(R0,T,0,2)(2,3:4)';
 R=[R1,R2]; eigenvalue=eig(R); mae=max(abs(eigenvalue));
 if(abs(c-0.7)<=1e-6) eig070=eigenvalue endif
 if(abs(c-0.65)<1e-6) eig065=eigenvalue; endif
 fpen_stab=[fpen_stab;c,p(1),mae];
endfor
plot(fpen_stab(:,1),fpen_stab(:,3), "k-", 'linewidth', 5); grid on;
xlabel("c"); ylabel("Largest |s_i|");
axis([0.64,0.7,0.6,1.4]);drawnow();
save fpen-stab.dat fpen_stab eig070 eig065
#print -deps2 -F:20 fpen-stab.eps

Code 35 “fpen-stabpf.m”

# "fpen-stabpf.m"
load fpen-bd.dat #fpen_bd; 分岐図
load fpen-stab.dat #fpen_stab eig070 eig065; 安定判別
eig070
eig065
### 安定性で分割 ###
i=1; mae=fpen_stab(i,3);
ps=[]; pu=[];
for i=1:rows(fpen_stab);
 mae=fpen_stab(i,3);
 if (mae<1)
  ps=[ps;fpen_stab(i,1:2)]; #安定
 else
  pu=[pu;fpen_stab(i,1:2)]; #不安定
 endif
endfor
hold on; pu2=pu(1:2:rows(pu),:); #表示を粗く
plot(ps(:,1),ps(:,2), "k-", 'linewidth', 2);
plot(pu2(:,1),pu2(:,2), "ko", 'markersize', 6);
plot(fpen_bd(:,1),fpen_bd(:,2), ".");
axis([0.61,0.68,-7.9,-6.5]);
xlabel("c"); ylabel("p_n"); drawnow();
#print -deps2 -F:20 fpen-stabpf.eps

Code 36 “fpen_lyap.cc”

Octave の開発環境で，Linux の場合，

> mkoctfile fpen_lyap.cc

として fpen_lyap.oct を作っておく．後述のプログラムで使う．

/* "fpen_lyap.cc"
 *
 * [...]$ mkoctfile fpen_lyap.cc
 *
 * として fpen_lyap.oct に変換する．
 */
#include <octave/oct.h>
#include <octave/LSODE.h>

/* 運動方程式
 * ※ octave でベクトルの第1成分は x(1) だが，ここでは x(0)
 * --- ここから --- */
struct _p { double c, k, P, omeg; } p;
ColumnVector xdot(const ColumnVector& x, double t) {
  ColumnVector dx(4);
  dx(0) = x(1); /* octave では dx(1)=x(2) */
  dx(1) = -p.c*x(1) -p.k*sin(x(0)) + p.P*cos(p.omeg*t);
  dx(2) = x(3); /* octave では dx(3)=x(4) */
  dx(3) = -p.k*cos(x(0))*x(2) -p.c*x(3) ;
  return dx;
}
/* --- ここまで --- */

double norm2(ColumnVector x) {
  double r=0.0;
  int dim=x.rows()/2;
  for (int i=dim; i<dim*2; i++) {
    r += x(i)*x(i);
  }
  return sqrt(r);
}
void normalize(ColumnVector &x, double r) {
  int dim=x.rows()/2;
  for (int i=dim; i<dim*2; i++) {
    x(i)/=r;
  }
}

DEFUN_DLD (fpen_lyap, args, , "Poincare map for a single pendlumn") {
  /* 引数の取得 */
  double get_global_double( char* s);
  ColumnVector x0 (args(0).vector_value ()); //初期値
  double T = (args(1).double_value ()); //周期
  double m = (args(2).int_value ()); //無視する過渡の点数
  double n = (args(3).int_value ()); //出力点数
  /* global 変数の取得 */
  p.c = get_global_double( "c" );
  p.k = get_global_double( "k" );
  p.P = get_global_double( "P" );
  p.omeg = get_global_double( "omeg" );
  /* 数値積分 */
  int dim = x0.rows(); //次元
  int Tn = 200; //分割数
  double t=0.0, dt = T/(double)Tn;
  ODEFunc odef(xdot);
  ColumnVector x=x0; //解
  x(2)=1.0; x(3)=0.0;
  double r = norm2(x);
  for ( int i=1; i<m; i++ ) { /* 過渡を無視する */
    normalize(x,r); /* 長さを1に */
    LSODE ls(x, 0.0, odef);
    for ( int ti=1; ti<=Tn; ti++ ) {
      t = dt*(double)ti;
      x = ls.do_integrate(t);
    }
    r = norm2(x);
  }
  Matrix pm(n,dim);
  ColumnVector lam(n-1);
  for ( int j=0; j<dim; j++ ) {
    pm(0,j)=x(j);
  }
  for ( int i=1; i<n; i++ ) {
    normalize(x,r); /* 長さを1に */
    LSODE ls(x, 0.0, odef);
    for ( int ti=1; ti<=Tn; ti++ ) {
      t = dt*(double)ti;
      x = ls.do_integrate(t);
    }
    r = norm2(x);
    lam(i-1)=log(r)/T;
    for ( int j=0; j<dim; j++ ) {
      pm(i,j)=x(j);
    }
  }
  octave_value_list retval;
  retval(1) = octave_value(lam);
  retval(0) = octave_value(pm);
  return retval;
}

double get_global_double( char* s ) {
  double retval = 0.123456789012345;
  octave_value tmp = get_global_value (s, true);
  if ( ! tmp.is_defined () || tmp.is_empty() )
    octave_stdout<< "Global variable \"" << s << "\" is empty";
  else
    retval = tmp.double_value ();
  return retval;
}

Code 37 “fpen-lyap.m”

※同じフォルダに，上で作った fpen_lyap.oct が必要．

# "fpen-lyap.m"
global c k P omeg;
c=0.3; k=1; P=2.7; omeg=1.0; T=2*pi/omeg;
n1=1000; #無視する過渡応答
n2=2000;  #データ数
x0=[0;0;1;1];
### カオスのリアプノフ指数 ###
[xx,lam]=fpen_lyap(x0,T,n1,n2);
lyap=mean(lam)
x=mod(xx(:,1).+pi,2*pi).-pi; #角度を-πからπに収める変換
plot(x,xx(:,2),"."); sleep(3);
plot(lam,'.;{/Symbol l_k};'); hold on;
plot(lyap*ones(size(lam)),'-;{/Symbol l};','linewidth',8);
xlabel("k"); ylabel("Lyapunov exponent");
drawnow(); hold off;
#print -deps2 -F:24 fpen-lyap.eps

Code 38 “fpen-lyapbd.m”

※同じフォルダに，上で作った fpen_lyap.oct が必要．

# "fpen-lyapbd.m"
global c k P omeg;
c=0.3; k=1; P=2.7; omeg=1.0; T=2*pi/omeg;
n1=800; #無視する過渡応答
n2=800;  #データ数
x0=[0;0;1;1];
### 周期倍加分岐のリアプノフ指数 ###
N=100; p=linspace(0.7,0.6,N+1);
hold on; lyapbd=[];
for i=1:N+1
 c=p(i); [xx,lam]=fpen_lyap(x0,T,n1,n2);
 lyap=mean(lam); graph=[p(i),lyap];
 plot(graph(1),graph(2),"."); drawnow();
 lyapbd=[lyapbd;graph];
 x0=xx(rows(xx),:)';
endfor
hold off; plot(lyapbd(:,1),lyapbd(:,2),"k-",'linewidth',5);
xlabel('c', 'fontsize', 32);
ylabel('Lyapunov exponent', 'fontsize', 32);
grid on; drawnow();
save fpen-lyapbd.dat lyapbd;
#print -deps2 -F:26 fpen-lyapbd.eps