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: 168
Registered: 4/4/09
Re: Automating an iterative procedure
Posted: Feb 25, 2013 3:39 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

It is working now, but how to implement convergence?
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;




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.