Module: AdaptiveIntegrator
Overview
This is the top-level driver, with logging, used by the more complex integrators in “Numerical Methods in C”, second edition. Code is adapted from C, of course.
Constant Summary collapse
- MAXSTP =
10000
- TINY =
1.0e-30
Instance Method Summary collapse
- #adaptive_integrate(ystart, x1, x2, eps, h1, hmin, derivs) ⇒ Object
- #integrate_ad_step(y, dydx, x, h, eps, yscal, derivs) ⇒ Object
Methods included from Integrator
#Integrator_initialize, #integrate, #set_max_samples
Instance Method Details
#adaptive_integrate(ystart, x1, x2, eps, h1, hmin, derivs) ⇒ Object
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 |
# File 'lib/integrator/adaptive.rb', line 28 def adaptive_integrate(ystart, x1, x2, eps, h1, hmin, derivs) x = x1 h1 = -h1 unless h1 >= 0 x2 - x1 >= 0 ? h = h1 : h = -h1 @nok = @nbad = @count = 0 y = ystart xsav = x - @dxsav * 2.0 if @kmax > 0 nstp = 1 while nstp <= MAXSTP dydx = derivs.call(x, y) #yscal = y.abs + (dydx * h).abs + Vector( [TINY] * y.length) yscal = y.collect2(dydx) { |yi, dyi| yi.abs + (dyi * h).abs + TINY } if @kmax > 0 and @count < @kmax - 1 and (x-xsav).abs > @dxsav.abs @count += 1 @xp[@count] = x @yp[@count] = y xsav = x end h = x2 - x if ((x + h - x2)*(x + h - x1) > 0.0) x, y, hdid, hnext = integrate_ad_step(y, dydx, x, h, eps, yscal, derivs) if hdid == h @nok += 1 else @nbad += 1 end if (x-x2)*(x2-x1) >= 0.0 if @kmax != 0 @count += 1 @xp[@count] = x @yp[@count] = y end return y end raise "Step size too small!" if hnext.abs <= hmin h = hnext nstp += 1 end raise "Too many steps in integration!" end |
#integrate_ad_step(y, dydx, x, h, eps, yscal, derivs) ⇒ Object
24 25 26 |
# File 'lib/integrator/adaptive.rb', line 24 def integrate_ad_step(y, dydx, x, h, eps, yscal, derivs) raise "Adaptive integrators must implement integrate_ad_step!" end |