ヤコビ法による固有値と固有ベクトルの算出(ウィウキンソンの公式使用)

function [v,e]=jacobiw(a)
n=size(a,1)
v=eye(a)
eps=0.000001*max(sum(abs(a),'c'))
kmax=400
for repeat=1:kmax
apqmax=0
for p=1:n
for q=1:n
if p<>q
apqmax=max(abs(a(p,q)),apqmax)
end
end
end
if apqmax<eps end
for p=1:n-1
for q=p+1:n
x=(a(p,p)+a(q,q))*0.5
y=(a(p,p)-a(q,q))*0.5
z=sqrt(y*y+a(p,q)*a(p,q))
if y<0 end
a(p,p)=x+z
a(q,q)=x-z
c=sqrt((1+y/z)*0.5)
s=a(p,q)*0.5/z/c
a(p,q)=0
a(q,p)=0
for j=1:n
if j<>p & j<>q
apj=a(p,j)
a(p,j)=apj*c+a(q,j)*s
a(q,j)=-apj*s+a(q,j)*c
a(j,p)=a(p,j)
a(j,q)=a(q,j)
end
end
for i=1:n
vip=v(i,p)
v(i,p)=vip*c+v(i,q)*s
v(i,q)=-vip*s+v(i,q)*c
end
end
end
e=diag(a)
end

LR 分解による固有値の算出

function e=lr(a)
clear crout
getf('crout.sci')
n=size(a,1)
m=100
for k=1:m
[l,r]=crout(a)
a=r*l
end
e=diag(a)

ヤコビ法による固有値と固有ベクトルの算出

function [v,e]=jacobi(a)
n=size(a,1)
v=eye(a)
eps=0.000001*max(sum(abs(a),'c'))
kmax=400
for repeat=1:kmax
apqmax=0
for p=1:n
for q=1:n
if p<>q
apqmax=max(abs(a(p,q)),apqmax)
end
end
end
if apqmax<eps end
for p=1:n-1
for q=p+1:n
app=a(p,p)
if app-a(q,q)<>0
theta=0.5*atan(2*a(p,q)/(app-a(q,q)))
else
theta=%pi*0.25
end
c=cos(theta)
s=sin(theta)
a(p,p)=app*c*c+2*a(p,q)*c*s+a(q,q)*s*s
a(q,q)=app*s*s-2*a(p,q)*c*s+a(q,q)*c*c
a(p,q)=0
a(q,p)=0
for j=1:n
if j<>p & j<>q
apj=a(p,j)
a(p,j)=apj*c+a(q,j)*s
a(q,j)=-apj*s+a(q,j)*c
a(j,p)=a(p,j)
a(j,q)=a(q,j)
end
end
for i=1:n
vip=v(i,p)
v(i,p)=vip*c+v(i,q)*s
v(i,q)=-vip*s+v(i,q)*c
end
end
end
e=diag(a)
end

以下のようにして出したものを使い
答え合わせをしてみましょう．

[v e]=spec(a)

ガウス・ザイデル法による連立1次方程式の解法

function x=seidel(a,b,x0,d)
//x0: 初期値
//d: 収束判定定数
d=d*d
n=length(x0)
x=x0
sa=1
while sa>d
for i=1:n
x(i)=b(i)
for j=1:n
if j<>i
x(i)=x(i)-a(i,j)*x(j)
end
end
x(i)=x(i)/a(i,i)
end
sa=0
x0=x-x0
sa=x0'*x0
x0=x
end

以下のようにして出したものを使い
答え合わせをしてみましょう．

a\\b

部分ピボット選択付きガウスの消去法を用いた連立1次方程式の解法

function x=gauss(a,b)
n=length(b)
for k=1:n-1
[res,ind]=sort(abs(a(k:n,k)))
saidai=ind(1)+k-1
if abs(a(saidai,k))<1e-7 end
if saidai<>k
kari=a(k,:)
a(k,:)=a(saidai,:)
a(saidai,:)=kari
kari=b(k)
b(k)=b(saidai)
b(saidai)=kari
end
a(k,k+1:n)=a(k,k+1:n)/a(k,k)
b(k)=b(k)/a(k,k)
a(k,k)=1
for i=k+1:n
a(i,k+1:n)=a(i,k+1:n)-a(i,k)*a(k,k+1:n)
b(i)=b(i)-a(i,k)*b(k)
a(i,k)=0
end
end
x(n)=b(n)/a(n,n)
for k=n-1:-1:1
x(k)=b(k)-a(k,k+1:n)*x(k+1:n)
end

クラウト法を用いて下三角行列と上三角行列に分解

function [l u]=crout(a)
n=size(a,1)
l=zeros(a)
u=zeros(a)
l(1,1)=1
u(1,:)=a(1,:)
l(:,1)=a(:,1)/u(1,1)
for i=2:n
l(i,i)=1
u(i,i)=a(i,i)-l(i,1:i-1)*u(1:i-1,i)
for j=i+1:n
l(j,i)=a(j,i)-l(j,1:i-1)*u(1:i-1,i)
l(j,i)=l(j,i)/u(i,i)
u(i,j)=a(i,j)-l(i,1:i-1)*u(1:i-1,j)
end
end

ロジックシステム

シミュレーション時間によって異なる信号を出力するシステムです．
図はシミュレーション時間を t とした場合，次のようになる例です．
1) 3<=t<=7 では，値 0 の一定信号を出力
2) 上記以外では，振幅 2，周波数 0.5 Hz の正弦波を出力

Block	Palette	Set	Value
Relational Op	Others	Operator	3
logical Op AND	Others
If then else	Branching
Gen\_Sin	Sources	Magnitute	2
		Frequency	\%pi

トリガ信号に同期

Block	Palette	Set	Value
Square generator	Sources
Edge trigger	Branching		1
Selector	Branching

移動平均

Block	Palette	Set	Value
Gen\_Sin	Sources	Magnitute	10
		Frequency	0.5*\%pi
Gen\_Sin	Sources	Frequency	40*\%pi
Clock	Sources	Period	0.01
sum	Linear	Number Of inputs	5
1/z	Linear	Inherit	1
Scope	Sinks	Ymin	-15 -15
		Ymax	15 15
		Refresh period	10

離散系

Block	Palette
S/H	Linear