线性代数学习二,用克拉默法则求解线性方程组

  这是书上1.4的内容,进度得加快

  1. #!/usr/bin/env lua
  2.  
  3. require 'determinant'
  4.  
  5. function Cramer(n)
  6.         assert(n>0)
  7.  
  8.         local cramer={}
  9.  
  10.         -- the coefficient determinant
  11.         cramer.D={}
  12.  
  13.         -- this D1,D2,D3...
  14.         cramer.Ds={}
  15.         for i=1,n do
  16.                 table.insert(cramer.Ds,{})
  17.         end
  18.  
  19.         -- this number of unknown numbers
  20.         cramer.n=n
  21.  
  22.         -- add an equation
  23.         -- for example:
  24.         -- x1+2x2-x3+3x4=2
  25.         -- use cramer:add(1,2,-1,3,2)
  26.         cramer.add=function(this,...)
  27.                 local toadd={...}
  28.                 assert(table.maxn(toadd)==this.n+1)
  29.  
  30.                 for i=1,this.n do
  31.                         table.insert(this.D,toadd[i])
  32.                         for j=1,this.n do
  33.                                 table.insert(this.Ds[j],toadd[i])
  34.                         end
  35.                 end
  36.  
  37.                 local last=toadd[this.n+1];
  38.                 local times=table.maxn(this.D)/this.n
  39.                 for i=1,this.n do
  40.                         this.Ds[i][(times-1)*this.n+i]=last
  41.                 end
  42.         end
  43.  
  44.         -- Use Cramer Rules to solve linear equations
  45.         -- When solution does not exist, nil is returned,
  46.         -- and solution table returned otherwise.
  47.         cramer.solve=function(this)
  48.                 assert(table.maxn(this.D)==this.n*this.n)
  49.  
  50.                 local D=Determinant(this.D):eval();
  51.                 if D==0 then
  52.                         return nil
  53.                 end
  54.  
  55.                 local solutions={}
  56.                 local Dn;
  57.                 for i=1,this.n do
  58.                         Dn=Determinant(this.Ds[i]):eval();
  59.                         table.insert(solutions,Dn/D)
  60.                 end
  61.  
  62.                 return solutions;
  63.         end
  64.  
  65.         return cramer
  66. end

下面是一个例子:

 

  1. cramer=Cramer(4)
  2. cramer:add(1,2,-1,3,2)
  3. cramer:add(2,-1,3,-2,7)
  4. cramer:add(0,3,-1,1,6)
  5. cramer:add(1,-1,1,4,-4)
  6.  
  7. solutions=cramer:solve()
  8. for _,v in ipairs(solutions) do
  9.         print(v)
  10. end

输出是1,3,2,-1,符合书上答案

Algorithm Comments(0) 2008年8月22日 11:00