Date: Mar 1, 2013 7:56 AM
Author: Milos Milenkovic
Subject: Random generator of initial matrix

>
Dear,
there is a procedure for iterative determining B0 and B1 matrices of 4x4 (It is given below). For the first iteration B0 is given and B1 have to be determined randomly, values have to be chosen from the interval [-0.1,0.1]. This first selection will greatly influence on the subsequent results of B0 and B1.
How to program a random generator which will generate values of matrix B1 in iteration 0 and use (save) those which will produce the suitable values for B0 and B1 in iteration 1. Suitable values are those within 0-1 interval.
Best,
Milos

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;