|
|
Re: Automating an iterative procedure
Posted:
Feb 25, 2013 3:39 AM
|
|
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;
|
|
