Finding the fixed points of a function

Question

I am trying to find the fixed point of a logistic distribution function and determine how the fixed point changes for different parameter values. The code looks like:

nfxp.reps <- 0
err <- 10
p <- seq(0, 1, by = 0.0001)
pold <- p
gamma <- 6
k <- 3
while (err > 1E-12) {
  nfxp.reps <- nfxp.reps + 1
  if (nfxp.reps > 999) { 
    stop("The number of NFXP reps needs to be increased. \n")
  } 
  pnew <- plogis(-k + gamma * pold) 
  err <- max(abs(pnew - pold))
  pold <- pnew
}

The above code works very well in the above parameter choices: gamma and k - find 3 fixed points, 2 stable and 1 unstable (where p=0.5). However if I change the above parameter non-proportionally, where the middle fixed point is either above or below 0.5, say for:

gamma<-7
k<-3

The loop is unable to locate the middle fixed point which is p=0.3225 (if gamma=7, k=3)

Answer 1

Fixed point iteration by construction cannot find the unstable equilibria in your setup since it is repelling. In other words, unless you start right at the unstable equilibria the nfxp algorithm will always move away from it.

An alternative approach is to use a root solving approach. Of course, there are no guarantees that all fixed points will be found. Here is a simple example:

library(rootSolve) # for the uniroot.all function
pfind<-function(k=3,gamma=7) 
{
pdiff <-function(p0) p0-plogis(-k + gamma * p0) 
uniroot.all(p.diff,c(0,1))
}
> fps= pfind()
> fps
[1] 0.08036917 0.32257992 0.97925817

We can verify this:

pseq =seq(0,1,length=100)
plot(x=pseq ,y= plogis(-k + gamma *pseq),type= 'l')
abline(0,1,col='grey')
points(matrix(rep(fps,each=2), ncol=2, byrow=TRUE),pch=19,col='red')

Hope this helps.

Answer 2

I rearrange your code in a new function.

p.fixed <- function(p0,tol = 1E-9,max.iter = 100,k=3,gamma=7,verbose=F){
  pold <- p0
  pnew <-  plogis(-k + gamma * pold) 
  iter <- 1
    while ((abs(pnew - pold) > tol) && (iter < max.iter)){
      pold <- pnew
      pnew <- plogis(-k + gamma * pold) 
      iter <- iter + 1
      if(verbose)
         cat("At iteration", iter, "value of p is:", pnew, "\n")
    }
    if (abs(pnew - pold) > tol) {
      cat("Algorithm failed to converge")
      return(NULL)
    }
    else {
      cat("Algorithm converged, in :" ,iter,"iterations \n")
      return(pnew)
    }
}

some tests:

p.fixed(0.2,k=3,gamma=7)
Algorithm converged, in : 30 iterations 
[1] 0.08035782
> p.fixed(0.2,k=5,gamma=5)
Algorithm converged, in : 7 iterations 
[1] 0.006927088
> p.fixed(0.2,k=5,gamma=5,verbose=T)
At iteration 2 value of p is: 0.007318032 
At iteration 3 value of p is: 0.006940548 
At iteration 4 value of p is: 0.006927551 
At iteration 5 value of p is: 0.006927104 
At iteration 6 value of p is: 0.006927089 
At iteration 7 value of p is: 0.006927088 
Algorithm converged, in : 7 iterations 
[1] 0.006927088

Answer 3

I don't really understand which distribution you are using; this is my standard code for the fixed point method which I always use and change if needed (you have to fill in your function f(x) in ftn;

fixed_point <- function(x0, eps = 1e-6, max_iter = 100){
  x.old <- x0
  x.new <- ftn(x.old)
  iter <- 1
  while((abs(x.new-x.old) > eps) && (iter < max_iter){
    x.old <- x.new
    x.new <- ftn(x.old)
    iter <- iter + 1
  }
 if (abs(x.new-x.old) > eps){
  cat("failed to converge\n")
  return(NULL)
 } else {
  return(x.new)
 }
}

Answer 4

Not sure what exactly you did wrong but I will give you my code which always works for finding the fixed point. The last function below can be used to compute function g, which is defined as g(x) = c*ftn(x) + x.

fixpt_own <- function(x0, tol = 1e-6, max.iter = 100) {
  xold <- x0
  xnew <- ftn_g(xold)
  iter <- 1
  cat("At iteration 1 value of x is:", xnew, "\n")
  while ((abs(xnew-xold) > tol) && (iter < max.iter)) {
    xold <- xnew;
    xnew <- ftn_g(xold);
    iter <- iter + 1
    cat("At iteration", iter, "value of x is:", xnew, "\n")
  }
  if (abs(xnew-xold) > tol) {
    cat("Algorithm failed to converge\n")
    return(NULL)
  } else {
    cat("Algorithm converged\n")
    return(xnew)
  }
}
fixpt_own(3,1e-6,150)


ftn_g <- function(x){
  c <- 4;
  g <- c*(((1+x)/x - log(x))/(1+x)^2) + x;
  return(g)
}