1-d Integrate, differentiate, solve ODE

spec/diff_spec.cr

@@ -0,0 +1,14 @@
+require "./spec_helper"
+
+describe GSL::Diff do
+  it "differentiate simple function" do
+    f = ->(x : Float64) { Math.log(x) }
+    result, eps = GSL.diff(f, 4.0)
+    result.should be_close 0.25, 1e-9
+  end
+
+  it "differentiate simple function in yield syntax" do
+    result, eps = GSL.diff(0.25) { |x| Math.log(x) }
+    result.should be_close 4.0, 1e-9
+  end
+end

spec/integration_spec.cr

@@ -0,0 +1,97 @@
+require "./spec_helper"
+
+describe GSL::Integration do
+  describe "qng algorithm" do
+    it "integrates simple function" do
+      f = ->(x : Float64) { Math.sin(x) }
+      result, eps, count = GSL::Integration.qng(f, 0.0, Math::PI/2)
+      result.should be_close 1.0, 1e-9
+      eps.should be < 1e-9
+      count.should be > 1
+    end
+
+    it "integrates closured function" do
+      a = 1.0
+      f = ->(x : Float64) { a*Math.sin(x) }
+      a = 8.0
+      result, _, _ = GSL::Integration.qng(f, 0.0, Math::PI/2)
+      result.should be_close 8.0, 1e-9
+    end
+
+    it "integrates with yielding syntax" do
+      a = 16.0
+      result, _, _ = GSL::Integration.qng(0.0, Math::PI) { |x| a*Math.sin(x) }
+      result.should be_close 32.0, 1e-9
+    end
+  end
+
+  describe "qag algorithm" do
+    it "integrates simple function" do
+      f = ->(x : Float64) { Math.sin(x) }
+      result, eps = GSL::Integration.qag(f, 0.0, Math::PI/2)
+      result.should be_close 1.0, 1e-9
+      result, eps = GSL::Integration.qag(f, 0.0, 10001*Math::PI)
+      result.should be_close 2.0, 1e-9
+    end
+
+    it "integrates function with singularity on a bound" do
+      f = ->(x : Float64) { 1.0/Math.sqrt(2 - x) }
+      result, eps = GSL::Integration.qags(f, 1.0, 2.0)
+      result.should be_close 2.0, 1e-9
+    end
+
+    it "integrates function with known singularity points" do
+      f = ->(x : Float64) { 1.0/Math.sqrt((2 - x).abs) }
+      result, eps = GSL::Integration.qagp(f, [1.0, 2.0, 3.0])
+      result.should be_close 4.0, 1e-9
+    end
+
+    it "integrates function on semi-infinite interval" do
+      # f = ->(x : Float64) { Math.exp(-x**2) }
+      f = ->(x : Float64) { 1.0 / (x**2 + 1) }
+      result, eps = GSL::Integration.qags(f, -Float64::INFINITY, Float64::INFINITY)
+      result.should be_close Math::PI, 1e-9
+
+      result, eps = GSL::Integration.qags(f, 0, Float64::INFINITY)
+      result.should be_close Math::PI/2, 1e-9
+
+      result, eps = GSL::Integration.qags(f, -Float64::INFINITY, 0)
+      result.should be_close Math::PI/2, 1e-9
+
+      expect_raises(ArgumentError) { GSL::Integration.qags(f, -Float64::INFINITY, -Float64::INFINITY) }
+      expect_raises(ArgumentError) { GSL::Integration.qags(f, Float64::INFINITY, -Float64::INFINITY) }
+      expect_raises(ArgumentError) { GSL::Integration.qags(f, Float64::INFINITY, 0.0) }
+      expect_raises(ArgumentError) { GSL::Integration.qags(f, 0.0, -Float64::INFINITY) }
+    end
+  end
+
+  describe "qawc algorithm" do
+    it "find Cauchy principal value for a singularity at finite point" do
+      f = ->(x : Float64) { 1.0 } # function is 1.0 / x
+      result, eps = GSL::Integration.qawc(f, -1, 1, 0)
+      result.should be_close(0.0, 1e-9)
+      f = ->(x : Float64) { Math.sqrt((2 - x).abs) } # function is 1.0 / Math.sqrt((2 - x).abs)
+      result, eps = GSL::Integration.qawc(f, 1, 4, 2)
+      result.should be_close(2*Math.sqrt(2) - 2, 1e-9)
+    end
+  end
+
+  describe "cquad algorithm" do
+    it "integrates simple function" do
+      f = ->(x : Float64) { Math.sin(x) }
+      result, eps, count = GSL::Integration.cquad(f, 0.0, Math::PI/2)
+      result.should be_close 1.0, 1e-9
+    end
+
+    it "integrates function with singularities" do
+      f = ->(x : Float64) { 1.0/Math.sqrt((2 - x).abs) }
+      result, eps, count = GSL::Integration.cquad(f, 1, 3)
+      result.should be_close 4.0, 1e-6
+    end
+  end
+
+  it "has high-level routine to integrate" do
+    r = GSL::Integration.integrate(0.0, Math::PI) { |x| Math.sin(x) }
+    r.should be_close 2.0, 1e-9
+  end
+end

spec/ode_spec.cr

@@ -0,0 +1,89 @@
+require "./spec_helper"
+
+class TestODE < GSL::ODE::System
+  def initialize(@k : Float64)
+    super(2)
+  end
+
+  def function(t, y, dydt)
+    dydt[0] = y[1]
+    dydt[1] = -@k * y[0]
+  end
+end
+
+class TestODEJacobian < GSL::ODE::JacobianSystem
+  def initialize(@mu : Float64)
+    super(2)
+  end
+
+  def function(t, y, dydt)
+    dydt[0] = y[1]
+    dydt[1] = -y[0] - @mu * y[1]*(y[0]*y[0] - 1)
+  end
+
+  def jacobian(t, y, dfdy, dfdt)
+    dfdy[0] = 0
+    dfdy[1] = 1.0
+    dfdy[2] = -2.0*@mu*y[0]*y[1] - 1.0
+    dfdy[3] = -@mu*(y[0]*y[0] - 1.0)
+    dfdt[0] = 0.0
+    dfdt[1] = 0.0
+  end
+end
+
+describe GSL::ODE do
+  describe "integrates system without jacobian" do
+    it "with natural step" do
+      sys = TestODE.new(1.0)
+      ode = GSL::ODE::Driver.new(sys, 1e-2, epsabs: 1e-9)
+      y0 = [1, 0]
+      y, t = ode.evolve(y0, 0, Math::PI/2)
+      y.last[0].should be_close(0.0, 1e-6)
+      y.last[1].should be_close(-1, 1e-6)
+      t.last.should be_close(Math::PI/2, 1e-6)
+      y, t = ode.evolve(y.last, Math::PI/2, Math::PI*2)
+      y.last[0].should be_close(1, 1e-6)
+      y.last[1].should be_close(0, 1e-6)
+      t.last.should be_close(Math::PI*2, 1e-6)
+    end
+
+    it "with given time values" do
+      sys = TestODE.new(1.0)
+      ode = GSL::ODE::Driver.new(sys, 1e-2, epsabs: 1e-9)
+      y0 = [0, 10]
+      results = ode.evolve(y0, [0, Math::PI/2, Math::PI, 3*Math::PI/2, Math::PI*2])
+      r = results.map { |y| y.map(&.round).to_a }
+      r.should eq [[0.0, 10.0], [10.0, 0.0], [0.0, -10.0], [-10.0, 0.0], [0.0, 10.0]]
+    end
+
+    it "with given time step yielding results" do
+      sys = TestODE.new(1.0)
+      ode = GSL::ODE::Driver.new(sys, 1e-2, epsabs: 1e-9)
+      y0 = [0, 10]
+
+      expected = [[0.0, 10.0], [10.0, 0.0], [0.0, -10.0], [-10.0, 0.0], [0.0, 10.0]]
+      i = 1
+      ode.evolve(y0, 0, Math::PI*2, Math::PI/2) do |y, t|
+        expected_y = expected[i]
+        expected_t = i*Math::PI/2
+        t.should be_close expected_t, 1e-6
+        y[0].should be_close expected_y[0], 1e-6
+        y[1].should be_close expected_y[1], 1e-6
+        i += 1
+      end
+    end
+  end
+  describe "integrates system with jacobian" do
+    it "with given time values" do
+      sys = TestODEJacobian.new(10.0)
+      ode = GSL::ODE::Driver.new(sys, 1e-2, epsabs: 1e-9, algorithm: GSL::ODE::Algorithm::MSBDF)
+      y0 = [1, 0]
+      y = ode.evolve(y0, [0, 9.5, 11.0, 19.0, 20.1])
+      y[0][0].should eq 1.0
+      y[1][0].should be_close -1.1, 1e-1
+      y[2][0].should be_close 2.0, 1e-2
+      y[3][0].should be_close 1.1, 1e-1
+      y[4][0].should be_close -2.0, 1e-2
+    end
+  end
+end

spec/poly_spec.cr

@@ -79,6 +79,17 @@ describe GSL do
       poly = GSL::Poly.new([1.0, 2.0, 1.0]) # x^2+2x+1
       poly.eval_derivs(0.1).should eq [1.21, 2.2, 2]
     end
+
+    it "#diff returns derivative of polynomial" do
+      poly = GSL::Poly.new([2.0, 1.0, 2.0, 1.0])         # x^3+2x^2+x+2
+      poly.diff.should eq GSL::Poly.new([1.0, 4.0, 3.0]) # 3*x^2+4x+1
+      poly.should eq GSL::Poly.new([2.0, 1.0, 2.0, 1.0])
+    end
+
+    it "#integrate returns integral of polynomial" do
+      poly = GSL::Poly.new([2.01, 1.0, 3.0, 1])                         # x^3+3x^2+x+2.01
+      poly.integrate.should eq GSL::Poly.new([0, 2.01, 0.5, 1.0, 0.25]) # x^4/4+x^3+x^2/2+2.01x
+    end
   end
 
   describe GSL::PolyDD do

spec/spec_helper.cr

@@ -1,2 +1,4 @@
 require "spec"
 require "../src/gsl"
+
+puts "GSL Version is: #{String.new(LibGSL.gsl_version)}"

src/gsl/base/function.cr

@@ -63,4 +63,26 @@ module GSL
   def self.wrap_function(&function : FunctionFDF)
     wrap_function(function)
   end
+
+  alias MultiRootFunction = (Tuple(Vector, Vector))
+  alias MultiRootFunctionFDF = (Vector -> Tuple(Float64, DenseMatrix))
+
+  # def self.wrap_function(function : MultiRootFunction)
+  #   result = uninitialized LibGSL::Gsl_multiroot_function
+  #   if function.closure?
+  #     box = Box.box(function)
+  #     # (LibC::Double, Void* -> LibC::Double)
+  #     result.f = ->(x : Gsl_vector*, data : Void*, f : Gsl_vector*) do
+  #       Box(MultiRootFunction).unbox(data).call(x, f)
+  #     end
+  #     result.params = box # no need to keep reference to `box` as `result` holds it.
+  #   else
+  #     result.params = function.pointer
+  #     result.function = ->(x : Float64, data : Void*) do
+  #       Function.new(data, Pointer(Void).null).call(x)
+  #     end
+  #   end
+  #   result
+  # end
+
 end

src/gsl/base/iterator.cr

@@ -1,7 +1,7 @@
 require "./matrix.cr"
 
 module GSL
-  class Vector
+  struct Vector
     def each(&block : Float64 -> _)
       self.to_a.each &block
       self

src/gsl/base/libgsl.cr

@@ -4135,7 +4135,7 @@ lib LibGSL
   fun gsl_linalg_complex_QR_unpack(qr : Gsl_matrix_complex*, tau : Gsl_vector_complex*, q : Gsl_matrix_complex*, r : Gsl_matrix_complex*) : LibC::Int
   fun gsl_linalg_complex_QR_unpack_r(qr : Gsl_matrix_complex*, t : Gsl_matrix_complex*, q : Gsl_matrix_complex*, r : Gsl_matrix_complex*) : LibC::Int
   fun gsl_linalg_QR_band_decomp_L2(m : LibC::SizeT, p : LibC::SizeT, q : LibC::SizeT, ab : Gsl_matrix*, tau : Gsl_vector*) : LibC::Int
-  fun gsl_linalg_QR_band_unpack_L2(p : LibC::SizeT, q : LibC::SizeT, qrb : Gsl_matrix*, tau : Gsl_vector*, q : Gsl_matrix*, r : Gsl_matrix*) : LibC::Int
+  fun gsl_linalg_QR_band_unpack_L2(p : LibC::SizeT, q : LibC::SizeT, qrb : Gsl_matrix*, tau : Gsl_vector*, qm : Gsl_matrix*, r : Gsl_matrix*) : LibC::Int
   fun gsl_linalg_QRPT_decomp(a : Gsl_matrix*, tau : Gsl_vector*, p : Gsl_permutation*, signum : LibC::Int*, norm : Gsl_vector*) : LibC::Int
   fun gsl_linalg_QRPT_decomp2(a : Gsl_matrix*, q : Gsl_matrix*, r : Gsl_matrix*, tau : Gsl_vector*, p : Gsl_permutation*, signum : LibC::Int*, norm : Gsl_vector*) : LibC::Int
   fun gsl_linalg_QRPT_solve(qr : Gsl_matrix*, tau : Gsl_vector*, p : Gsl_permutation*, b : Gsl_vector*, x : Gsl_vector*) : LibC::Int
@@ -5019,7 +5019,7 @@ lib LibGSL
     params : Void*
   end
 
-  fun gsl_multiroot_fdjacobian(f : Gsl_multiroot_function*, x : Gsl_vector*, f : Gsl_vector*, epsrel : LibC::Double, jacobian : Gsl_matrix*) : LibC::Int
+  fun gsl_multiroot_fdjacobian(func : Gsl_multiroot_function*, x : Gsl_vector*, f : Gsl_vector*, epsrel : LibC::Double, jacobian : Gsl_matrix*) : LibC::Int
   fun gsl_multiroot_fsolver_alloc(t : Gsl_multiroot_fsolver_type*, n : LibC::SizeT) : Gsl_multiroot_fsolver*
 
   struct Gsl_multiroot_fsolver_type
@@ -5492,7 +5492,7 @@ lib LibGSL
   fun gsl_ran_lognormal_pdf(x : LibC::Double, zeta : LibC::Double, sigma : LibC::Double) : LibC::Double
   fun gsl_ran_logarithmic(r : Gsl_rng*, p : LibC::Double) : LibC::UInt
   fun gsl_ran_logarithmic_pdf(k : LibC::UInt, p : LibC::Double) : LibC::Double
-  fun gsl_ran_multinomial(r : Gsl_rng*, k : LibC::SizeT, n : LibC::UInt, p : LibC::Double*, n : LibC::UInt*)
+  fun gsl_ran_multinomial(r : Gsl_rng*, k : LibC::SizeT, n : LibC::UInt, p : LibC::Double*, nn : LibC::UInt*)
   fun gsl_ran_multinomial_pdf(k : LibC::SizeT, p : LibC::Double*, n : LibC::UInt*) : LibC::Double
   fun gsl_ran_multinomial_lnpdf(k : LibC::SizeT, p : LibC::Double*, n : LibC::UInt*) : LibC::Double
   fun gsl_ran_negative_binomial(r : Gsl_rng*, p : LibC::Double, n : LibC::Double) : LibC::UInt

src/gsl/base/vector.cr

@@ -1,5 +1,5 @@
 module GSL
-  class Vector
+  struct Vector
     getter pointer
 
     def initialize(@size : Int32)

src/gsl/maths/analysis/diff.cr

@@ -0,0 +1,27 @@
+module GSL
+  enum Diff::Direction
+    Central
+    Forward
+    Backward
+  end
+
+  def self.diff(function : GSL::Function, x : Float64, step : Float64 = 0.01, dir : Diff::Direction = Diff::Direction::Central)
+    f = wrap_function(function)
+    result = uninitialized Float64
+    abserr = uninitialized Float64
+    case dir
+    in .central?
+      code = LibGSL.gsl_deriv_central(pointerof(f), x, step, pointerof(result), pointerof(abserr))
+    in .forward?
+      code = LibGSL.gsl_deriv_forward(pointerof(f), x, step, pointerof(result), pointerof(abserr))
+    in .backward?
+      code = LibGSL.gsl_deriv_backward(pointerof(f), x, step, pointerof(result), pointerof(abserr))
+    end
+    code = LibGSL::Code.new(code)
+    return result, abserr
+  end
+
+  def self.diff(x : Float64, step : Float64 = 0.01, dir : Diff::Direction = Diff::Direction::Central, &f : GSL::Function)
+    diff(f, x, step, dir)
+  end
+end