Search All of the Math Forum:

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

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: Random generator of initial matrix
Replies: 0

 Milos Milenkovic Posts: 189 Registered: 4/4/09
Random generator of initial matrix
Posted: Mar 1, 2013 7:56 AM

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