sidi.R

This implement two polynomial methods to accelerate general fixed point iterations. See http://www.stat.ucla.edu/~deleeuw/janspubs/2008/reports/deleeuw_R_08i.pdf
language: R
license: GPL 2

Code for Snippet:

                
# A general procedure to accelerate fixed point iterations by using either MPE or RRE
# proedures. See Sidi, J. Comp. Sci, 3, 2012, 92-101
 
 
sidi<-function(f,xinit,r,itmax=100,eps=1e-10,expol=mpe,verbose=FALSE,noisy=FALSE,stable=TRUE) {
xold<-xinit; itel<-1; xx<-matrix(0,r+2,length(xinit)); oaps<-1; oerr<-1; aps<-Inf; s<-1/(r+1)
repeat{
    xx[1,]<-xold
    for (i in 2:(r+2)) xx[i,]<-f(xx[i-1,])
    if (r == 0) xnew<-xx[2,]
        else xnew<-expol(xx,r)
    xstb<-f(xnew)
    naps<-norm(xold-xnew); crat<-naps/oaps
    nerr<-norm(xnew-xstb); erat<-nerr/oerr
    if (stable) xnew<-xstb
    if (verbose) {
        cat("Iteration: ",
            formatC(itel,digits=6,width=6),
            "   Change:   ",
            formatC(naps,digits=10,width=15,format="f"),
            "   Error:   ",
            formatC(nerr,digits=10,width=15,format="f"),
            "  CRatio:   ",
            formatC(crat,digits=10,width=15,format="f"),
            "  ERatio:   ",
            formatC(erat,digits=10,width=15,format="f"),
           "\n")
        }
    if (noisy) {
       cat("Solution:  ",
            format(xnew,digits=10,width=15,format="f"),
           "\n")
        }
    if ((itel == itmax) || (nerr < eps)) break()
    itel<-itel+1
    xold<-xnew
    oaps<-naps
    oerr<-nerr
    }
return(list(itel=itel,aps=naps,err=nerr,crat=crat,erat=erat,x=xnew))
}
 
mpe<-function(xx,r) {
u<-diff(xx)
if (r == 1) uu<-as.matrix(u[1,])
    else uu<-t(u[1:r,])
cc<-c(lsfit(uu,-u[r+1,],intercept=FALSE)$coef,1) 
cc/sum(cc)   
return(drop(cc%*%xx[-(r+2),]))
}
 
rre<-function(xx,r) {
u<-t(diff(xx)); v<-rowMeans(u)
ac<-lsfit(u[,-(r+1)]-u[,r+1],-u[,r+1],intercept=FALSE)$coef
cc<-c(ac,1-sum(ac))
return(drop(cc%*%xx[-(r+2),]))  
}
 
comments powered by Disqus

Info

Source Site:

Jan

Tags: Convergence Acceleration Fixed Points

Link to this snippet:


Download to Code Collector

To use the direct link to your snippet on CodeCollector.net either copy the html from the above section or drag the Download to Code Collector to where you would like to use it.

More Info:

Times Viewed: 603
Date Added: 2013-03-06 23:16:43
Last Modified: 2013-04-17 22:12:19

Web Analytics