Thread Subject:

Gram Schmidt orthogonalization

Hi, I've been trying to implement Gram Schmidt orthogonalization for calculating lyapunov exponents. However, I find that I slowly lose orthogonality over about 1000 time steps; I think it may be due to normalization after each orthogonalization step. Below is my code:

 %Gram-Schmidt orthogonalization and normalization
  orthcol=Jacobianmatrix;
  orthcol(:,1)=orthcol(:,1)/sqrt(transpose(orthcol(:,1))*(orthcol(:,1)));

    for jj=2:3 %for 2nd and 3rd vectors
        for kk=1:jj-1
           orthcol(:,jj)=orthcol(:,jj)-dot(orthcol(:,jj),orthcol(:,kk))...
               /sqrt(transpose(orthcol(:,kk))*(orthcol(:,kk)))*orthcol(:,kk);
        end
         orthcol(:,jj) = orthcol(:,jj)/sqrt(transpose(orthcol(:,jj))*(orthcol(:,jj))); %normalize
        
    end

Any help would be really appreciated. Thanks!

My previous program implemented the classical Gram-Schmidt and that didn't work either.

"Andrew Greene" wrote in message <ik7ppc$ie2$1@fred.mathworks.com>...
> Hi, I've been trying to implement Gram Schmidt orthogonalization for calculating lyapunov exponents. However, I find that I slowly lose orthogonality over about 1000 time steps; I think it may be due to normalization after each orthogonalization step. Below is my code:

Yes, it is well known Gram Schmidt is slightly affected by numerical error. Use QR() that carry out the orthogonalisation in a stable way.

Bruno

hi bruno,

Thanks for the reply. I was also using QR but I had wanted to code something of my own; is this not possible due to the errors that creep in? I had thought it would be possible because I have been given code in C that does not suffer from the loss of orthogonality.

Also, do you know how I would extract expansion/contraction information for my Lyapunov exponents as QR automatically orthonormalizes?

best regards,

AG.

"Andrew Greene" wrote in message <ik8673$svb$1@fred.mathworks.com>...
> hi bruno,
>
> Thanks for the reply. I was also using QR but I had wanted to code something of my own; is this not possible due to the errors that creep in? I had thought it would be possible because I have been given code in C that does not suffer from the loss of orthogonality.

There is another algorithm so called "modified Gram Schmidt" which is slightly more stable than the original one.
 
>
> Also, do you know how I would extract expansion/contraction information for my Lyapunov exponents as QR automatically orthonormalizes?

Sorry, I don't. I only know Lyapunov exponents by name and by the theorem, and never work closely with them.

Bruno