|
|
Re: Automating an iterative procedure
Posted:
Feb 23, 2013 11:22 AM
|
|
Solution is based on for iterator and symplifying the code by making of functions for certain parts.
Best,
M
"Milos Milenkovic" <m.milenkovic@mathworks.com> wrote in message <kg7clr$166$1@newscl01ah.mathworks.com>...
> %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!!
|
|
