View original (French)

Cholesky Method Implementation in MATLAB

serhanus Posted messages 3 Status Member -  
serhanus Posted messages 3 Status Member -
Hello everyone,

First of all, I remind you that the Cholesky method consists of solving a classical linear system Ax=b, by decomposing the matrix A into A=B*Bt (with B a lower triangular matrix, and Bt its transpose).

There is a function in MATLAB called "chol" that returns Bt; I wrote a script that solves the system using this function.

But the problem is that I am asked to implement a function that does what "chol" does!

1 answer

  1.
    rlo73 Posted messages 3123 Status Member 669
     
    It's called an exercise, isn't it? We're not going to write the algorithm for you, after all, it seems to me that's not CCM's policy!
    3
    1. serhanus Posted messages 3 Status Member
       
      yes it is,
      but I'm not asking you to solve the problem, the problem is that it's the first publication, and it was posted before completion...

      what I want is just an indication

      here is my attempt:

      function x=cholesky2(A,b)
      %*****************
      [n nn]=size(A);
      %*****************
      for i=1:n
      for j=1:n
      B(i,j)=0;
      end
      end

      %*****************
      for j=1:n-1
      som1=0;

      for k=1:n
      som1=som1+B(j,k)^2;
      end

      B(j,j)=sqrt(A(j,j)-som1);

      n=j+1;
      som2=0;
      for k=1:j-1
      som2=som2+B(j,k)*B(j+1,k);
      end
      B(j+1,j)=A(j+1,j)-som2;
      end
      %B(j+1,j)=A(j+1,j)-som2;
      B
      %*****************
      Bt=B';
      %******************
      for k=1:n
      B(k,n+1)=b(k);
      end

      %******************
      y(1)=B(1,n+1)/B(1,1);
      for i=1:n;
      c=0;

      for j=1:i-1;


      c=c+B(i,j)*y(j);


      end
      y(i)=(B(i,n+1)-c)/B(i,i);
      end
      %******************

      for k=1:n
      Bt(k,n+1)=y(k);
      end


      %******************
      x(n)=Bt(n,n+1)/Bt(n,n);
      for i=n-1:-1:1;
      c=0;

      for j=i+1:n;


      c=c+Bt(i,j)*x(j);


      end
      x(i)=(Bt(i,n+1)-c)/Bt(i,i);
      end
      x=x';
      %******************
      B
      Bt
      0