0

So I am attempting to implement an algorithm that uses Newton's method for finding the root of some given function, and I am having what I imagine is decent success except at one small but very frustrating point.

When I evaluate my function for small initial values of x_0 (i.e. less than or close to 1, [-1.1,1.1]), my function seems to work fine, but when I attempt to use it for larger values (even something like 1.2 or higher) i get this error message:

Error in while (abs(x_1 - x_0) > 10^(-1 * k) && iter < max_iter) { : 
  missing value where TRUE/FALSE needed

Here is my code for this problem:

g_x <- function(x){
  return(x/(1+x^2)^(1/2))
}

dg_x <- function(x){
  # derivative of the given function of g
  return(1/(x^2+1)^(3/2))
}

newton <- function(x_0,k){
  # x_0 is initial guess from which you want to being zeroing in on the root
  # k is the parameter of the tolerance = 10^(-k)
  if(between(g_x(x_0),-10^(-1*k)/2,10^(-1*k)/2)){
    return(c(x_0,0))
  }
  
  iter <- 1
  x_1 <- x_0 - (g_x(x_0)/dg_x(x_0))
  
  while(abs(x_1-x_0) > 10^(-1*k)){
    x_0 <- x_1
    x_1 <- x_0 - (g_x(x_0)/dg_x(x_0))
    iter <- iter + 1
  }
  return(c((x_1+x_0)/2, iter))
}

I have no idea what i need to do to fix this. I have been going through this line by line in the console, and the computer is capable of evaluating the expression abs(x_1-x_0) > 10^(-k) as TRUE for the first couple iterations. I am not sure why it becomes a problem at all.

It is not like the computer reches some super small/large value of x_1 that it generalizes as zero or infinity because this problem presents itself on the very first iteration of the algorithm.

In the problem we are given that: x_0 = 2 and g(x) = x/sqrt(x^2+1) so dg(x) = 1/(x^2+1)^3/2

therefore the next point to calculate is x_1 = -8

so the computer is evaluating abs(-8-2) > 10^-k. Why does this break

jay.sf
  • 60,139
  • 8
  • 53
  • 110
Chance Law
  • 1
  • Just put a print inside the while loop and you'll see that you are overflowing. – nicola Dec 02 '22 at 12:52
  • Take also a look at https://en.wikipedia.org/wiki/Newton%27s_method#Proof_of_quadratic_convergence_for_Newton's_iterative_method. It looks that over 1 you are not in the convergence range. This is due to the fact that the term proportional to the second derivative in the Taylor expansion is bigger. If you start at 1, you stay there forever. So, if math is not an opinion, you are not going to fix this. – nicola Dec 02 '22 at 13:08
  • On a basic level, this is a duplicate of https://stackoverflow.com/q/7355187/3358272. Troubleshooting may be easier if your error message and your code agreed, but they don't: the error cites `(abs(x_1 - x_0) > 10^(-1 * k) && iter < max_iter)` but the code only has `(abs(x_1-x_0) > 10^(-1*k))`, and while the `&& ...` might not be a player, is there anything else we're missing. – r2evans Dec 02 '22 at 13:22
  • Suggestion: add _before_ your `while` statement `test <- (abs(x_1-x_0) > 10^(-1*k)); if (length(test) != 1L || anyNA(test)) browser();`; add the same thing to the end of the `while` block, and then change the while itself to `while (test)`. This won't fix your problem, but when you run it, you will be placed into a debugger when the condition causing the error is met. Once there, check the values of all `x_1`, `x_0`, and `k` to see which of them is not what you expect. From there, you'll need to trace back further in the function to see where it went wrong. – r2evans Dec 02 '22 at 13:25

1 Answers1

0

Before doing a non-linear optimization problem, I will typically look at the function that I am seeking to solve with a graph. Since your function only depends on one variable it is straightforward to do this.

xs = seq(-10,10,length=10000)
plot(xs,g_x(xs))

yields

enter image description here

So the 'root' is at 0 (the function has a value of 0 at 0), but the derivative is "very spiky" which may mean that the Newton approach will repeatedly over shot the solution repeatedly.

I modified your function to put in an upper and lower bound for x, and it works provided that you start reasonably close to the solution:

newton <- function(x_0,k,lower=-10,upper=10){
  # x_0 is initial guess from which you want to being zeroing in on the root
  # k is the parameter of the tolerance = 10^(-k)
  if(between(g_x(x_0),-10^(-1*k)/2,10^(-1*k)/2)){
    return(c(x_0,0))
  }
  
  iter <- 1
  x_1 <- x_0 - (g_x(x_0)/dg_x(x_0))
  
  while( (abs(x_1-x_0) > 10^(-1*k)) & between(x_0,lower,upper)){
    x_0 <- x_1
    x_1 <- x_0 - (g_x(x_0)/dg_x(x_0))
    iter <- iter + 1
  }
  return(c((x_1+x_0)/2, iter))
}

#x_0=10
#k=10
newton(.3,10)

newton(5,10)

Starting at 0.3 works (solves for -3.812802e-15 ), and starting at 5 does not (exits at 976500). Using a coarse grid search to choose a good starting value for your function would be a reasonable approach.

This shows what happens if you start at say 2:

x_0=2  
iter <- 1
x_1 <- x_0 - (g_x(x_0)/dg_x(x_0))

for (inter in 1:10){
  x_0 <- x_1
  x_1 <- x_0 - (g_x(x_0)/dg_x(x_0))
  iter <- iter + 1
  print(x_1,iter)
}

returns:

[1] 512
[1] -1.34e+08
[1] 2.418e+24
[1] -1.4135e+73
[1] 2.82401e+219
[1] NaN
[1] NaN
[1] NaN
[1] NaN
[1] NaN
DashdotdotDashdotdot
  • 449
  • 1
  • 9