|
|
Re: Automating an iterative procedure
Posted:
Feb 23, 2013 3:04 PM
|
|
Dear,
is this possible:
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];
for i=1:100
B1(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];
B0(1)=[-1 0 0 0; 0 -1 0 0 ; 0 0 -1 0; 0 0 0 -1];
Ke[i]=wn1(B0(i),B1(i),K0,K1); %wn d1 d2 up are approapriate functions
[Delta1(i), Error1(i)]=d1(B1(i),K1,Ke(i));
[Delta0(i), Error0(i)]=d2(B0(i),B1(i),K0,Delta1(i),Ke(i));
[B0(i),B1(i)]= up(B0(i),B1(i), Delta0(i), Delta1(i))
end
or i,j indexing?
"Milos Milenkovic" <m.milenkovic@mathworks.com> wrote in message <kgaqbg$853$1@newscl01ah.mathworks.com>...
> 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!!
|
|
