newtonplot.jl

using Plots, ForwardDiff, StaticArrays, Distributions, Test

g(x) = sin.(x)

function newton(g, x0, n=10)
    out = zeros(length(x0), n)
    x = x0
    for i in 1:n
        out[:, i] = x
        dg = ForwardDiff.jacobian(g, x)
        x = x - dg \ g(x)
    end
    return out
end

function quasinewton(g, x0, n=10)
    out = zeros(length(x0), n)
    x = x0
    dg = ForwardDiff.jacobian(g, x)
    for i in 1:n
        out[:, i] = x
        x = x - dg \ g(x)
    end
    return out
end

function newtonplot(g, x0; n=10, op=newton, title="newton's method", xstar=0)
  n = 10

  xs = op(g, [x0], n)
  ys = g.(xs)
  xs = xs[:]
  ys = ys[:]

  p = plot(sin, range(-3.0, 3.0, length=100), xlim=(-pi, pi), legend=false, title=title, foreground_color_border=:transparent, foreground_color_axis=:transparent)

  for i in 1:n-1
      plot!(p, Shape([ (xs[i], 0), (xs[i], ys[i]) ]), linecolor=:orange)
      plot!(p, Shape([ (xs[i], ys[i]), (xs[i+1], 0) ]))
  end
  plot!(p, [xs[1]], [0.], marker=true, markerstrokewidth=0)

  if abs(xstar - xs[n]) < .01
    plot!(p, [xstar], [0.], marker=true, markercolor=RGB(.3,.9,0.), markerstrokewidth=0, markersize=5.)
  else
    plot!(p, [xstar], [0.], marker=true,
  xs = op(g, [x0], n)
  ys = g.(xs)
  xs = xs[:]
  ys = ys[:]

  p = plot(sin, range(-3.0, 3.0, length=100), xlim=(-pi, pi), legend=false, title=title, foreground_color_border=:transparent, foreground_color_axis=:transparent)

  for i in 1:n-1
      plot!(p, Shape([ (xs[i], 0), (xs[i], ys[i]) ]), linecolor=:orange)
      plot!(p, Shape([ (xs[i], ys[i]), (xs[i+1], 0) ]))
  end
  plot!(p, [xs[1]], [0.], marker=true, markerstrokewidth=0)

  if abs(xstar - xs[n]) < .01
    plot!(p, [xstar], [0.], marker=true, markercolor=RGB(.3,.9,0.), markerstrokewidth=0, markersize=5.)
  else
    plot!(p, [xstar], [0.], marker=true, markercolor=:red, markerstrokewidth=0, markersize=5.)
  end

  p
end

anim = @animate for y in range(0, 2pi, length=180)
  x = cos(y) * 1.3
  plot(newtonplot(g, x, op=newton), newtonplot(g, x, op=quasinewton, title="quasinewton method"), layout=(2,1), dpi=200)
end
        
gif(anim, "newton.gif")