Scicos の C block を使う (LSI C-86)

LSI C-86
を
C:Program Filesscilab-4.1.2lcc
に入れます．

角度制御ステップ応答 (ロボコンマガジン46) Scicos

最急降下法を用いて非線形関数の最小値を求める

f(x1,x2)=(x1-2)^2+(x1-2*x2)^2
x1, x2 の初期値はどちらも 0
ステップ幅は 0.1

alpha=0.1;
m=100;
x=zeros(2,1);
for i=1:m
dx(1)=2*(x(1)-2)+2*(x(1)-2*x(2));
dx(2)=-4*(x(1)-2*x(2));
if norm(dx)<%eps
break
end
x=x-alpha*dx;
printf("x1=%f, x2=%fn",x(1),x(2));
end
printf("x1=%f, x2=%f,minimum=%fn",x(1),x(2),...
(x(1)-2)^2+(x(1)-2*x(2))^2);

ニュートン法による非線形方程式の解法

f1=2*x1^3-x2^3-9=0
f2=x1^3-2*x2^3-4=0

を初期値 (0.5,1) で解いています．

x=[0.5;1];
for i=1:20
f(1)=2*x(1)*x(1)*x(1)-x(2)*x(2)*x(2)-9;
f(2)=x(1)*x(1)*x(1)-2*x(2)*x(2)*x(2)-4;
j(1,1)=6*x(1)*x(1);
j(1,2)=-3*x(2)*x(2);
j(2,1)=3*x(1)*x(1);
j(2,2)=-6*x(2)*x(2);
dx=j\f;
if norm(dx)<%eps
break
end
x=x-dx;
printf("x1=%f, x2=%f\n",x(1),x(2));
end

角度制御ステップ応答 (ロボコンマガジン46)

clear
deg2rad=%pi/180;
rad2deg=180/%pi;
tt=0:0.001:2;
thref=ones(tt)*60*deg2rad;
Kp=0.2;
Kd=0.025;
J=0.001;
s=poly(0,'s');
G1=syslin('c',1,J*s);
G2=syslin('c',1,s);
Dcont=syslin('c',Kd,1);
Pcont=syslin('c',Kp,1);
Vclosed=G1/.Dcont;
Xclosed=(Pcont*Vclosed*G2)/.syslin('c',1,1)
[th,XX]=csim(thref,tt,Xclosed);
plot(tt,th*rad2deg,tt,thref*rad2deg);
legend('th[deg]','thref[deg]');
xgrid;
a=gca();
a.data_bounds=[0 0;1.5 80];
xtitle('','time[s]','th[deg]');

ルンゲ・クッタ法による微分方程式の解法

以下は2階斉次線形微分方程式

d^2x/dt^2-dx/dt-2x-6=0

を

dx_1/dt=x_2
dx_2/dt=x_2+2x_1+6

のような連立1次微分方程式として解く場合の例です．

function [t,x]=runge(h)
n=0.2/h
t(1)=0
x(1,:)=[0 0]
for i=1:n
k1=h*f(t(i),x(i,:))
k2=h*f(t(i)+h/2,x(i,:)+k1/2)
k3=h*f(t(i)+h/2,x(i,:)+k2/2)
k4=h*f(t(i)+h,x(i,:)+k3)
x(i+1,:)=x(i,:)+(k1+2*k2+2*k3+k4)/6
t(i+1)=t(i)+h
end
endfunction

function dif=f(t,x)
dif(1,1)=x(2)
dif(1,2)=x(2)+2*x(1)+6
endfunction

オイラー法を用いた微分方程式の解法

以下は

dx/dt=-x

を解く場合の例です．

function [t,x]=euler(h)
n=2/h
t(1)=0
x(1)=1
for i=1:n
x(i+1)=x(i)+h*f(t(i),x(i))
t(i+1)=t(i)+h
end
endfunction

function dif=f(t,x)
dif=-x
endfunction

ラグランジュ多項式を用いた補間

function p=lagrange(x,y)
n=length(x)
s=poly(0,'s')
p=0
for i=1:n
q=1
for j=1:n
if j<>i
q=q*(s-x(j))/(x(i)-x(j))
end
end
p=p+q*y(i)
end

線形補間

function [xval,yval]=senkeihokan(x,y)
n=length(x)
div=10
xval=[]
for i=1:n-1
h=(x(i+1)-x(i))/div
xval=[xval x(i):h:x(i)+h*(div-1)]
yval((i-1)*div+1:i*div)=y(i)+(y(i+1)-y(i))* ...
(xval((i-1)*div+1:i*div)-x(i))/(x(i+1)-x(i))
end
xval=[xval x(n)]
yval=[yval y(n)]

セル配列と構造体

セル配列
a=makecell([2 2],'Keio',[1 2;3 4],ones(2,2),rand(3,3))
a(1,2)
b=cell2mat(a(1,2))
b(2,1)
c=cell(2,2)
c(1,1)=makecell(1,[1 2;3 4])
c(2,2)=makecell(1,[4 3;2 1])

構造体
student.name='korondemo'
student.score=[12 34 56]
student=struct('name','korondemo','score',[12 34 56])
student(2).name='yamada'
student(2).score=[78 90 100]
student.score
student(2).score(2)