Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

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


Math Forum » Discussions » Software » comp.soft-sys.matlab

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

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Milos Milenkovic

Posts: 149
Registered: 4/4/09
Re: Automating an iterative procedure
Posted: Feb 23, 2013 11:22 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

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!!




Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.