The general form of the algorithm is here:
$${b_{n} = b_{n-1}- {{f(x)}\over{f'(x)}}}$$
$$end:f(x) < err$$
It is simply an arbitrary starting value less the function of the value over the derivative of the value.The idea is that the result of the first iteration feeds into the subsequent iteration, and the process terminates when f(x) < err threshold.
We start with parameter b0 and an error threshold which remains constant.
We want to find the x intercept of the following function
$$f(x)=-x^3+2x^2-e^x$$
We will need its derivative
$${{df(x)\over{dx}} =-3x^2+4x-e^x}$$
The form of the algorithm will look as follows:
$${b_{1} = b_{0}- {{-x^3+2x^2-e^x}\over{-3x^2+4x-e^x}}}$$
We arbitrarily set b0 as 10 and the error threshold to .001
I will print the values out instead of computing the algorithm with MathJax.
b0: -6.47058884010279 fx: 354.649469421832 fxdx: -151.489463487155
b0: -4.12950543233229 fx: 104.509232053146 fxdx: -67.6925579111767
b0: -2.58562475655285 fx: 30.5816407170309 fxdx: -30.474214161426
b0: -1.58209959409794 fx: 8.76059234819741 fxdx: -14.0430588411127
b0: -0.95826169330173 fx: 2.33291052253774 fxdx: -6.97140224182995
b0: -0.623621624975194 fx: 0.484337035554275 fxdx: -4.19719802152472
b0: -0.508226298600212 fx: 0.046298105996643 fxdx: -3.4093487299717
b0: -0.494646547843232 fx: 0.000591761687825487 fxdx: -3.32239821278166
-0.4946465
The value of x where f(x) = 0 is ~ -.495. This can be seen in the above graph where the blue line intercepts the x axis.
R Script
library(ggplot2) #Function f1_x <- function(x) { return (-x^3+2*(x^2)-exp(x)) } #Derivative f1_dxdy <- function(x) { return (-3*(x^2)+(4*x)-exp(x)) } #Load data dfFX <- data.frame(x=seq(-2.5,2.5,.01),f1_x=f1_x(seq(-2.5,2.5,.01)), f1_dxdy=f1_dxdy(seq(-2.5,2.5,.01))) #plot the data ggplot(data=dfFX) + geom_line(aes(x=x, y=f1_x),size=1,color="blue") + ylim(c(-2.5,2.5)) + xlim(c(-2.5,2.5)) + geom_vline(xintercept=0) + geom_hline(yintercept=0) #The Newton Method - caps N NewtonM <- function(b0,err,fx=f1_x,fxdx=f1_dxdy) { b0 <- b0-(fx(b0)/fxdx(b0)) print(paste("b0:",b0,"fx:",fx(b0),"fxdx:",fxdx(b0))) if (abs(f1_x(b0)) <= err) { return(b0) } if (abs(f1_x(b0)) >= err) { return(newtonM(b0,err)) } }
No comments:
Post a Comment