Posted to tcl by egavilan at Sun May 17 22:05:22 GMT 2009view raw

  1. package require math::linearalgebra
  2.  
  3. set x {0 1 2 3 4 5 6 7 8 9 10}
  4. set y {1 6 17 34 57 86 121 162 209 262 321}
  5.  
  6. proc build.matrix {xvec degree} {
  7. set sums [llength $xvec]
  8. for {set i 1} {$i <= 2*$degree} {incr i} {
  9. set sum 0
  10. foreach x $xvec {
  11. set sum [expr {$sum + pow($x,$i)}]
  12. }
  13. lappend sums $sum
  14. }
  15.  
  16. set order [expr {$degree + 1}]
  17. set A [math::linearalgebra::mkMatrix $order $order 0]
  18. for {set i 0} {$i <= $degree} {incr i} {
  19. set A [math::linearalgebra::setrow A $i [lrange $sums $i $i+$degree]]
  20. }
  21. return $A
  22. }
  23.  
  24. proc build.vector {xvec yvec degree} {
  25. set sums [list]
  26. for {set i 0} {$i <= $degree} {incr i} {
  27. set sum 0
  28. foreach x $xvec y $yvec {
  29. set sum [expr {$sum + $y * pow($x,$i)}]
  30. }
  31. lappend sums $sum
  32. }
  33.  
  34. set x [math::linearalgebra::mkVector [expr {$degree + 1}] 0]
  35. for {set i 0} {$i <= $degree} {incr i} {
  36. set x [math::linearalgebra::setelem x $i [lindex $sums $i]]
  37. }
  38. return $x
  39. }
  40.  
  41.  
  42. # build the system A.x=b
  43. set degree 2
  44. set A [build.matrix $x $degree]
  45. set b [build.vector $x $y $degree]
  46. # solve it
  47. set coeffs [math::linearalgebra::solveGauss $A $b]
  48. # show results
  49. puts $coeffs