109  方阵A的Crout分解                                                                          
 1  功能  对方阵A作LU分解,其中L为下三角阵, U为单位上三角阵.
 2  数学方法简介  对奇异方阵A通常是消去主对角线以下的元素,使A变为下三角阵,记未变换的A为A(0)   .
    第一步,  若aij(0) 0, 令uij =aij(0) /a11(0)   , j=2,3,…..,n .构造初等方阵    
                  1    u12    u13 .  .  .u1n                     
              U1=        1    
                               1
                                     .
                                       1
则                               
                     1    -u12    -u13  . . . –u1n
             U1(-1) =          1    
                                    1    
                                            1
用u1(-1)  右乘A,
                         a11(0)    0      0      0
             AU1(-1)  =   a21(0)   a22(1)   a23(1) … a2n(1)
                                 …   …                              
                         an1(0)   an2(1)   an3(1) … a nn(1)
   第二步,若a22(1) 0, 令u2j =a2j(1)  /a22(1)  , j=3,4,……,n.  构造初等方阵
                     1     0       0  …   0
          u2(-1)=             1      -u23   -u2n
                                    1    
                                           1
用u2(-1) 右乘A(1),则
                              a11(0)    0      0  …   0
      A(1)U2(-1)=A(0) U1(-1)U2(-1)=   a21(0)   a22(0)    0  …    0
                              a31(0)  a32(1)   a33(2)   a3n(2)
                                  …   …   …
                              an1(0)  an2(1)   an3(2)   ann(2)
若a33(2)0,可按上述步骤继续往下作,akk(k-1)  0,k=1,2,……,n,作完n-1步,A U1(-1) U2(-1)…Un-1(-1)=A(n-1) =L
则L为下三角阵,令U(-1) =U1(-1)U2(-1)…Un-1(-1) ,它为单位上三角阵,则
            A=A(n-1) U=LU
这里L为下三角阵,U为单位上三角阵.
    为了省去一些工作单元,用紧凑格式编写程序,即将L存放在A的左下方(含对角线),将单位上三角阵U存放在A的右上方(不含对角线元,对角线元全为1不必存放).
上述变换过程的计算方法步骤为,对I=1,2,…,n, 有
      lji =aji –∑ljk uki  ,  j=i,i+1,…,n
      uij=(aij –∑ljk ukj)/ lii ,    j=i +1,i+2,…,n
每步先算i列的元素,再算i行的元素.
3   使用说明
输入参数
N       整变量,方阵A的阶数.
A(N,N)   N*N个元素的二维实数组,按行存放矩阵A.
输出参数
W       标志,W=0分解完毕,W=1分解进行不下去.
L(N,N)   下三角阵L;
U(N,N)   单位上三角阵U.
10   ‘*************************************’
20   ‘*   109   方阵的A的CROUT分解    *’
30   ‘*************************'************’
40   INPUT "N,E="; N, E
50   DIM A(N, N)
60   PRINT TAB(3); "EXAMPLE"; TAB(8); "JUZHEN A": PRINT
70   FOR I = 1 TO N
80     FOR J = 1 TO N
90     READ A(I, J): PRINT USING; "####.#"; A(I, J);
100      NEXT J
110      PRINT
120    NEXT I
130   PRINT: GOSUB 500
140   PRINT TAB(5); "W="; W: PRINT
150   IF (W = 1) THEN GOTO 380
160   PRINT TAB(15); "JUZHEN L"
170   PRINT
180   FOR I = 1 TO N
190   FOR J = 1 TO I
200   PRINT USING; "##.#####"; A(I, J);: PRINT " ";
210   NEXT J
220   PRINT
230   NEXT I
240   FOR I = 1 TO N
250   A(I, I) = 1
260   NEXT I
270   PRINT
280  PRINT TAB(15); "JUZHEN U"
290  PRINT
300  FOR I = 1 TO N
310  FOR J = I TO N
320  PRINT TAB(10 * (J - 1)); USING; "##.#####"; A(I, J);
330  NEXT J
340  PRINT
350  NEXT I
360  DATA 1, 2, 3, 4, 1, 4, 9, 16, 1, 8, 27, 64, 1, 16, 81, 256
370  DATA
380  END
500 'SON'
510 IF ABS(A(1,1))〈=E THEN GOTO 750
520 FOR J=2 TO N
530 A(1,J)=A(1,J)/A(1,1)
540 NEXT J
550 P=0
560 FOR K=2 TO N
570  FOR I=K TO N
580  FOR 1TO K-1
590  P=P+A(I,R)*A(R,K)
600    NEXT R
610  A(I,K)= A(I,K)-P
620  P=0
630  NEXT I
640   P=0 
650  FOR J=K +1TO N
660  FOR R=1TO K-1
670    P=P+A(K,R)*A(R,J)
680  NEXT R
690 IF ABS(A(K,K))〈=E THEN GOTO 750
700  A(K,J)=(A(K,J)-P)/A(K,K)
710  P=0
720  NEXT J
730  NEXT K
740  W=0:RETURN
750        W=1:RETURN
====================================================================
要求用C语言或者pascal做出来,请大家务必帮忙。 

解决方案 »

  1.   

    已经有人作出一部分了,请大家把后面的继续完成即可,谢谢了
    program JuZhen;{$APPTYPE CONSOLE}uses
      SysUtils;const
                            //360  DATA 1, 2, 3, 4, 1, 4, 9, 16, 1, 8, 27, 64, 1, 16, 81, 256
                            //370  DATA
      OriData:array[0..15] of Integer=(1, 2, 3, 4, 1, 4, 9, 16, 1, 8, 27, 64, 1, 16, 81, 256);function Tab(N:Integer):String;
    var
      i:Integer;
    begin
      Result:='';
      for i:=1 to N do
        Result:=Result+' ';
    end;var
      N,E,W:Integer;
      A:array of array of Integer;
      i,j:Integer;
    begin
      Writeln('N,E=');      //40   INPUT "N,E="; N, E
      Readln(N);
      Readln(E);
      SetLength(A,N);       //50   DIM A(N, N)
      for i:=0 to N-1 do
        SetLength(A[i],N);
                            //60   PRINT TAB(3); "EXAMPLE"; TAB(8); "JUZHEN A": PRINT
      Writeln(Tab(3)+'EXAMPLE'+Tab(8)+'JUZHEN A');
      for i:=0 to N-1 do    //70   FOR I = 1 TO N
      begin
        for j:=0 to N-1 do  //80     FOR J = 1 TO N
        begin               //90     READ A(I, J): PRINT USING; "####.#"; A(I, J);
          A[i,j]:=OriData[i*N+j];
          Write(Format('%4.1f',[A[i,j]*1.0]));
        end;                //100      NEXT J
        Writeln;            //110      PRINT
      end;                  //120    NEXT I                        //130   PRINT: GOSUB 500
      //GOSUB 500 ?? ******
                            //140   PRINT TAB(5); "W="; W: PRINT
      Writeln(Tab(5)+'W='+IntToStr(W));
      if W=1 then           //150   IF (W = 1) THEN GOTO 380
        exit;
                            //160   PRINT TAB(15); "JUZHEN L"
                            //170   PRINT
      Writeln(Tab(15)+'JUZHEN L');
      for i:=0 to N-1 do    //180   FOR I = 1 TO N
      begin
        for j:=0 to i do    //190   FOR J = 1 TO I
                            //200   PRINT USING; "##.#####"; A(I, J);: PRINT " ";
          Write(Format('%2.5f ',[A[i,j]*1.0]));
                            //210   NEXT J
        Writeln;            //220   PRINT
      end;                  //230   NEXT I
      for i:=0 to N-1 do    //240   FOR I = 1 TO N
        A[i,i]:=1;          //250   A(I, I) = 1
                            //260   NEXT I
      Writeln;              //270   PRINT
                            //280  PRINT TAB(15); "JUZHEN U"
                            //290  PRINT
      Writeln('               JUZHEN U');
      for i:=0 to N-1 do    //300  FOR I = 1 TO N
      begin
        for j:=i to N-1 do  //310  FOR J = I TO N
                            //320  PRINT TAB(10 * (J - 1)); USING; "##.#####"; A(I, J);
          Writeln(Tab(10*j)+Format('%2.5f',[A[i,j]*1.0]));
                            //330  NEXT J
        Writeln;            //340  PRINT
      end;                  //350  NEXT Iend.                    //380  END
      

  2.   

    就是把那段Basic翻译成C语言或者pascal就行。
      

  3.   

    //LU分解
    //A为n x n,index为n+1个元素的整形数组
    //结束时A为复合的LU分解阵,即有上三角阵的全部和下三角阵的除去主对角线部分,下三
    //角的主对角线为1,index返回行交换的信息,d为正负1,表示行交换为奇数还是偶数,
    //用于求行列式时确定符号
    //部分选主元,最后要用index把行交换回来
    procedure LU_Decompose(var A:TMatrix;var index:IVector;var d:Extended);
    var
      i,imax,j,k,n:Integer;
      big,dum,sum,temp:Extended;
      vv:array of Extended;
    begin
      imax:=0;
      if A.nRow<>A.nCol then
        raise Exception.Create('LU_Decompose:Invalid Matrix A');
      n:=A.nRow;
      SetLength(vv,n+1);
      SetLength(index,n+1);
      //for i:=1 to n do index[i]:=i;
      d:=1;
      for i:=1 to n do
      begin
        big:=0;
        for j:=1 to n do
        begin
          temp:=abs(A[i,j]);
          if temp>big then
            big:=temp;
        end;
        if big=0 then
          raise Exception.Create('LU_Decompose: Singular Matrix');
        vv[i]:=1.0/big;
      end;
      for j:=1 to n do
      begin
        for i:=1 to j-1 do
        begin
          sum:=A[i,j];
          for k:=1 to i-1 do
            sum:=sum-A[i,k]*A[k,j];
          a[i,j]:=sum;
        end;
        big:=0;
        for i:=j to n do
        begin
          sum:=a[i,j];
          for k:=1 to j-1 do
            sum:=sum-A[i,k]*A[k,j];
          a[i,j]:=sum;
          dum:=vv[i]*abs(sum);
          if dum>=big then
          begin
            big:=dum;
            imax:=i;
          end;
        end;
        if j<>imax then
        begin
          for k:=1 to n do
          begin
            dum:=a[imax,k];
            A[imax,k]:=A[j,k];
            A[j,k]:=dum;
          end;
          d:=-d;
          vv[imax]:=vv[j];
        end;
        index[j]:=imax;
        if A[j,j]=0 then
          A[j,j]:=TINY;
        if j<>n then
        begin
          dum:=1/A[j,j];
          for i:=j+1 to n do
            A[i,j]:=A[i,j]*dum;
        end;
      end;
    end;
      

  4.   

    补充一下,TMatrix是一个1基准的矩阵类,换成2维数组也可以
    nRow,nCol返回TMatrix的行数和列数
    type
      IVector=array of Integer;
    const
      TINY=1e-20;
      

  5.   

    不选主元的话。。
    //A为输入矩阵,A_L返回下三角,A_U返回上三角
    //因为不选主元,所以有的矩阵会失败(抛出DivideByZero的异常)
    procedure LU_Decompose(const A:TMatrix;var A_L,A_U:TMatrix);
    var
      n,i,j,k,l:Integer;
      sum:Extended;
    begin
      sum:=0;
      if A.nRow<>A.nCol then
        raise Exception.Create('LU_Decompose:Invalid Matrix A');
      n:=A.nRow;
      A_L.nRow:=n;
      A_L.nCol:=n;
      A_U.nRow:=n;
      A_U.nCol:=n;
      for i:=1 to n do
      begin
        A_L[i,i]:=1;
        for j:=i+1 to n do
          A_L[i,j]:=0;
        for j:=1 to i-1 do
          A_U[i,j]:=0;
      end;  
      for j:=1 to n do
      begin
        for i:=1 to j do
        begin
          sum:=0;
          if i<>1 then
            for k:=1 to i-1 do
              sum:=sum+A_L[i,k]*A_U[k,j];
          A_U[i,j]:=A[i,j]-sum;
        end;
        for i:=j+1 to n do
        begin
          sum:=0;
          for k:=1 to i-1 do
            sum:=sum+A_L[i,k]*A_U[k,j];
          A_L[i,j]:=(A[i,j]-sum)/A_U[j,j];
        end;
      end;
    end;