root_finding.fs

let bisect f l r =
  let rec bisect l r i =
    seq {
      let xl, xr = f l, f r
      let sl, sr = int (sign xl), int (sign xr)
      if sl <> sr then
        printfn "[%i] f(%f) = %f, f(%f) = %f" i l xl r xr
        yield (l, r)
        let c = (l + r) / 2.0
        let x = f c
        let s = int (sign x)
        if sl = s then
          yield! bisect c r (i+1)
        else
          yield! bisect l c (i+1)
      else
        while true do yield (l, r)
    }
  bisect l r 0

let falseposition f l r =
  let rec fp l r i =
    seq {
      let xl, xr = f l, f r
      let sl, sr = int (sign xl), int (sign xr)
      if sl <> sr then
        printfn "[%i] f(%f) = %f, f(%f) = %f" i l xl r xr
        yield (l, r)
        let c =
          l + abs (l - r) * abs xl / (abs xl + abs xr)
        let x = f c
        let s = int (sign x)
        if sl = s then
          yield! fp c r (i+1)
        else
          yield! fp l c (i+1)
      else
        while true do yield (l, r)
    }
  fp l r 0

let secant f (l: float) (r: float) =
  let rec s l r i =
    seq {
      let fl, fr = f l, f r
      printfn "[%i] f(%f) = %f, f(%f) = %f" i l fl r fr
      yield (l, r)
      let df = (fr - fl) / (r - l)
      let next = r - fr/df
      yield! s r next (i+1)
    }
  s l r 0

let newton f df x =
  let rec n x i =
    seq {
      let fx = f x
      let dfx = df x
      printfn "[%i] f(%f) = %f, f'(%f) = %f" i x fx x dfx
      yield x
      let next = x - fx/dfx
      yield! n next (i+1)
    }
  n x 0

let f x = x ** 3.0 - 3.0 * x + 1.0

// falseposition f 1.0 2.0 |> Seq.take 100 |> Seq.last |> printfn "%A"
// secant f 1.0 2.0 |> Seq.take 10 |> Seq.last |> printfn "%A"

let g x = x ** 3.0 + 3.0 * x + cos x
let dg x = 3.0 * (x ** 2.0) - sin x + 3.0

newton g dg (-2.0) |> Seq.take 6 |> Seq.iteri (printfn "%i: %A")