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
Automating an iterative procedure
Posted: Feb 22, 2013 4:10 AM

%Dear,
%I need help with automating the following iterative procedure:
%i=1,...,N; N is number of iterations, it depends on the convergence of
%B0[i] and B1[i] but let assume for now that it is given.
%Input parameters:
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];
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];
%for i=1:
B0[i=1]=[-1 0 0 0; 0 -1 0 0 ; 0 0 -1 0; 0 0 0 -1];
B1[i=1]=[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];
%Unknown parameters:
%for i=1:N =>
Ke0[i]; Ke1[i]; Ke[i];
B0[i+1]=B0[i]+Delta0[i];
B1[i+1]=B1[i]+Delta1[i];
%Iteration i:
condB1[i] = cond(B1[i]);
condB0[i] = cond(B0[i]);
condB1[i] = 1.4662e+017;
condB0[i] = 1;
Ke1[i] = pinv(B1[i])*K1/B0[i]';
Ke0[i]=dlyap(B0[i],K0,[],B1[i]);
Ke[i]=0.5*[Ke1[i]+Ke0[i]);
A1[i]=B1[i]*Ke[i]*B1[i]';
A2[i]=Ke[i]*B1[i]';
A3[i]=B1[i]*Ke[i];
A4[i]=K1-A2[i];
n=length(A2[i]);
P1[i] = zeros(n);
P1[i](1:end)=1:n^2;
Q1[i] = kron(speye(n),A3[i]);
Q1[i] = Q1[i](:,P1[i]');
M1[i] = kron(A2[i]',speye(n)) + Q1[i];
Delta1[i]= M1[i] \ A4[i](:);
Delta1[i] = reshape(Delta1[i],[n n]);
norm( Delta1[i]*A2[i]+A3[i]*Delta1[i]'-A4[i]) / norm(A4[i]);
A5[i]=B0[i]*Ke[i]*B0[i]';
A6[i]=B0[i]*Ke[i];
A7[i]=Ke[i]*B0[i]';
A8[i]=K0-A5[i]-A1[i]-A3[i]*Delta1[i]'-Delta1[i]*A2[i];
n=length(A7[i]);
P0[i] = zeros(n);
P0[i](1:end)=1:n^2;
Q0[i]= kron(speye(n),A6[i]);
Q0[i] = Q0[i](:,P0[i]');
M0[i] = kron(A7[i]',speye(n)) + Q0[i];
Delta0[i] = M0[i] \ A8[i](:);
Delta0[i] = reshape(Delta0[i],[n n]);
norm( Delta0[i]*A7[i]+A6[i]*Delta0[i]'-A8[i]) / norm(A8[i]);
%for the next iteration
B0[i+1]=B0[i]+Delta0[i];
B1[i+1]=B1[i]+Delta1[i];
and the iteration again.
%Thanks!!

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