Vendor tcllib linear algebra

1 files changed, 86 insertions, 0 deletions

diff --git a/vendor/math/calculus.testscript b/vendor/math/calculus.testscript
new file mode 100644
index 00000000..81dd091a
--- /dev/null
+++ b/vendor/math/calculus.testscript
@@ -0,0 +1,86 @@
+# calculus.test --
+#    Test cases for the Calculus package
+#
+source calculus.tcl
+
+#
+# Simple test functions - exact result predictable!
+#
+proc const_func { x } {
+   return 1
+}
+proc linear_func { x } {
+   return $x
+}
+proc downward_linear { x } {
+   return [expr {100.0-$x}]
+}
+proc downward_linear { x } {
+   return [expr {100.0-$x}]
+}
+
+#
+# Test the Integral proc
+#
+puts "[::Calculus::Integral 0 100 100 const_func] - expected: 100"
+puts "[::Calculus::Integral 0 100 100 linear_func] - expected: 5000"
+puts "[::Calculus::Integral 0 100 100 downward_linear] - expected: 5000"
+puts "[::Calculus::Integral 0 100 100 downward_linear] - expected: 5000"
+puts "[::Calculus::IntegralExpr 0 100 100 {100.0-$x}] - expected: 5000"
+
+proc const_func2d { x y } {
+   return 1
+}
+proc linear_func2d { x y } {
+   return $x
+}
+puts "[::Calculus::Integral2D { 0 100 10 } { 0 50 1 } const_func2d] - \
+ expected 5000"
+puts "[::Calculus::Integral2D { 0 100 1  } { 0 50 1 } const_func2d] - \
+ expected 5000"
+puts "[::Calculus::Integral2D { 0 100 10 } { 0 50 1 } linear_func2d] - \
+ expected 250000"
+
+# xvec should one long!
+proc const_func { t xvec } { return 1.0 }
+
+# xvec should be two long!
+proc dampened_oscillator { t xvec } {
+   set x  [lindex $xvec 0]
+   set x1 [lindex $xvec 1]
+   return [list $x1 [expr {-$x1-$x}]]
+}
+
+foreach method {EulerStep HeunStep} {
+   puts "Method: $method"
+
+   set xvec   0.0
+   set t      0.0
+   set tstep  1.0
+   for { set i 0 } { $i < 10 } { incr i } {
+      set result [::Calculus::$method $t $tstep $xvec const_func]
+      puts "Result ($t): $result"
+      set t      [expr {$t+$tstep}]
+      set xvec   $result
+   }
+
+   set xvec   { 1.0 0.0 }
+   set t      0.0
+   set tstep  0.1
+   for { set i 0 } { $i < 20 } { incr i } {
+      set result [::Calculus::$method $t $tstep $xvec dampened_oscillator]
+      puts "Result ($t): $result"
+      set t      [expr {$t+$tstep}]
+      set xvec   $result
+   }
+}
+
+#
+# Boundary value problems:
+# use simple functions
+#
+proc coeffs { x } { return {1.0 0.0 0.0} }
+proc forces { x } { return 0.0 }
+
+puts [::Calculus::BoundaryValueSecondOrder coeffs forces {0.0 1.0} {100.0 0.0} 10]
+puts [::Calculus::BoundaryValueSecondOrder coeffs forces {0.0 0.0} {100.0 1.0} 10]