Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Topic: Automating an iterative procedure
Replies: 9   Last Post: Feb 25, 2013 3:39 AM

 Messages: [ Previous | Next ]
 Milos Milenkovic Posts: 189 Registered: 4/4/09
Re: Automating an iterative procedure
Posted: Feb 24, 2013 2:02 PM

Sorry for toppost, I give the all procedure bellow. I think that I implemented all your suggestions. How to make n iterations? Maybe to avoid for loops by vectorization? In that case, without any indexing of B0 and B1, how to connect iteration i with these parameters in code?
Best,
M
K1=[-0.2582 -18.596 0 -0.2582; -34.737 -1300.188 0 -34.737; 0 0 0 0; -0.2582 -18.596 0 -0.2582];
K0=[2.337 48.551 0 2.337; 48.551 3072.487 0 48.551; 0 0 0 0; 2.337 48.551 0 2.337];
B1=[0.01 0.02 0.07 0.03; 0.02 0.05 0.06 0.02; 0.07 0.09 0.01 0.03; 0.01 0.03 0.02 0.04];
B0=[-1 0 0 0; 0 -1 0 0 ; 0 0 -1 0; 0 0 0 -1];

for i=1:100
Ke=wn1(B0,B1,K0,K1);
[Delta1, Error1]=d1(B1,K1,Ke);
[Delta0, Error0]=d2(B0,B1,K0,Delta1,Ke);
[B0,B1]= up(B0,B1, Delta0, Delta1)
end

wn1 ->
function [Ke]=wn1(B0,B1,K0,K1)
condB1 = cond(B1); condB0 = cond(B0);
condB1 = 1.4662e+017; condB0 = 1;
X = pinv(B1)*K1/B0';
%Y=A*X*B
E = norm(K1-B1*X*B0'); F = norm(B1);
G = norm(B1*X*B0'); relerr = E/F;
Y = dlyap(B1,K0,[],B0); Ke=0.5*(X+Y);

d1 ->
function [Delta1,Error1]=d1(B1,K1,Ke)
A2=B1*Ke*B1'; A21=Ke*B1';
A22=B1*Ke; K11=K1-A2;
n=length(A21); P1 = zeros(n);
P1(1:end)=1:n^2;
Q1 = kron(speye(n),A22);
Q1 = Q1(:,P1');
M = kron(A21',speye(n)) + Q1;
Delta1= M \ K11(:);
Delta1 = reshape(Delta1,[n n]);
Error1 = norm( Delta1*A21+A22*Delta1'-K11) / norm(K11)

d2->
function [Delta0, Error0]=d2(B0,B1, K0,Delta1, Ke)
A1=B1*Ke*B1'; A2=Ke*B1'; A3=B1*Ke;
A5=B0*Ke*B0'; A6=B0*Ke; A7=Ke*B0';
A8=K0-A5-A1-A3*Delta1'-Delta1*A2;
n=length(A7); P0 = zeros(n);
P0(1:end)=1:n^2;
Q0= kron(speye(n),A6);
Q0 = Q0(:,P0');
M0 = kron(A7',speye(n)) + Q0;
Delta0 = M0 \ A8(:);
Delta0 = reshape(Delta0,[n n]);
Error0= norm( Delta0*A7+A6*Delta0'-A8) / norm(A8)

up ->
function [B0,B1]= up(B0,B1, Delta0, Delta1)
B0=B0+Delta0;
B1=B1+Delta1;

Date Subject Author
2/22/13 Milos Milenkovic
2/23/13 Milos Milenkovic
2/23/13 Milos Milenkovic
2/23/13 dpb
2/23/13 Milos Milenkovic
2/23/13 dpb
2/24/13 Milos Milenkovic
2/24/13 dpb
2/24/13 Milos Milenkovic
2/25/13 Milos Milenkovic