Method: Interval#newton

Defined in:
lib/interval.rb

#newton(f, fp, opts = {}) ⇒ Object

Solve a non-linear equation using the Newton-Raphson method, finding all solutions in a in interval.

For instance, to solve

x**3 == x

we rewrite it first in the form

x**3 - x == 0

The left-hand side of this equation has derivative

3* x**2 - 1

To find all solutions in [-100,100] we call

Interval[-100,100].newton(proc {|x| x**3 - x}, proc {|x| 3* x**2 - 1})
                         # => Interval[[-1.0], [-0.0], [1.0]]

a sharp result showing the three roots -1, 0 and +1.



396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
 
# File 'lib/interval.rb', line 396

def newton(f, fp, opts = {})
  effectiveOpts = NewtonOptions.dup.update(opts)
  verbose = effectiveOpts[:verbose]
  puts "Starting on #{self}" if verbose
  step = proc {|w,ww| (w - f.call(Interval[w]) / fp.call(ww) ) & ww}
  self.map{|xx|
    effectiveOpts[:maxIterations].times{
      previous = xx;
      xx = step.call(xx.midpoint,xx)
      if previous == xx
        if xx.sharp?
          puts "Sharp fixed point #{xx}" if verbose
          xx = Interval[xx.extrema.select {|x| f.call(Interval[x]).include?(0)}]
          break
        end
        nonminimal_extrema = xx.extrema.reject {|x| f.call(Interval[x]).include?(0)}
        if nonminimal_extrema == [] then
          puts "Unsharp fixed point #{xx}" if verbose
          break
        end
        if nonminimal_extrema.each {|x|
            yy = step.call(x,xx)
            if yy != xx
              xx = yy;
              break
            end
          }
          xx = Interval[
            if nonminimal_extrema.include?(xx.inf)
              xx.inf + xx.inf.ulp
            else
              xx.inf
            end,
            if nonminimal_extrema.include?(xx.sup)
              xx.sup - xx.sup.ulp
            else
              xx.sup
            end]
        end
      end
      if xx.empty?
        puts "No solution" if verbose
        break
      elsif ! xx.simple?
        puts "Branch splitting" if verbose
        xx = xx.newton(f,fp,opts)
        break
      end
    } && verbose && puts("Failed convergence #{effectiveOpts[:maxIterations]}")
    xx
  }. inject{|a,x| a |=x}
end