2010年8月27日金曜日

2ブロックシステムのプログラム

// 2ブロック

function Xdot=f(t,X)

P=X(1)
A=X(2)
Z=X(3)
D=X(4)


//------------------------------------------------------

Xmat1=[
0, A*P/(1+P), 0, 0;
0, 0, 0.01*Z*0.1, 0.01*Z*0.9 + 0.01*A;
0, 0, 0, 0.1*Z;
0.01*D, 0, 0, 0
]

dX1=[
-sum( Xmat1(1,:) ) + sum( Xmat1(:,1) );
-sum( Xmat1(2,:) ) + sum( Xmat1(:,2) );
-sum( Xmat1(3,:) ) + sum( Xmat1(:,3) );
-sum( Xmat1(4,:) ) + sum( Xmat1(:,4) );

]

//------------------------------------------------------

P=X(5)
A=X(6)
Z=X(7)
D=X(8)

//------------------------------------------------------

Xmat2=[
0, A*P/(1+P), 0, 0;
0, 0, 0.01*Z*0.1, 0.01*Z*0.9 + 0.01*A;
0, 0, 0, 0.1*Z;
0.01*D, 0, 0, 0
]

dX2=[
-sum( Xmat2(1,:) ) + sum( Xmat2(:,1) );
-sum( Xmat2(2,:) ) + sum( Xmat2(:,2) );
-sum( Xmat2(3,:) ) + sum( Xmat2(:,3) );
-sum( Xmat2(4,:) ) + sum( Xmat2(:,4) );
]

Xdot=[dX1 * 0.9 + dX2 * 0.1;
dX1 * 0.1 + dX2 * 0.9]

endfunction

// 初期値

X1s=[2; 1; 0.01; 0]
X2s=[2; 3; 0.01; 1]

Xs=[X1s; X2s]

//時間

t0=0; t1=50; dt=0.01;
t= t0 : dt : t1;

//微分方程式solver
x=ode(Xs,t0,t,f)

// グラフの表示
clf()

scf(0)
plot2d(t,x(1,:),2)
plot2d(t,x(2,:),3)
plot2d(t,x(3,:),4)
plot2d(t,x(4,:),5)

scf(1)
plot2d(t,x(5,:),6)
plot2d(t,x(6,:),7)
plot2d(t,x(7,:),8)
plot2d(t,x(8,:),9)

0 件のコメント:

コメントを投稿