R中牛顿迭代法的编写

求解方程的数值解，已经有了第一部计算，现在的问题是 迭代计算 ，因为迭代计算我不太懂，所以程序也不会编
f(x)=0
x(k+1)=x(k)-{J(x(k))}^(-1)*f(x(k)), k=0,1,2….
其中J和f已经计算出来了，这个Newtons迭代函数编写如下？
Newtons<-function (funs, y, ep=1e-5, it_max=100){
index<-0; k<-1
while (k<=it_max){
y1 <- y; obj <- funs(y);
y<- y-solve(obj$J, obj$f);
。。。

其中要迭代的函数：
funs<-function(x){
。。。
f<- c(。。,。。。) #这是f(.)
J<- matrix(,..,….,…,
nrow=2, byrow=T)
# Jacobi矩阵 j(.)
return(
list(f=f, J=J))
}
要迭代部分函数是正确的，运行问题在于：
错误于drop(.Call("La_dgesv", a, as.matrix(b), tol, PACKAGE = "base")) :
  Lapack例行程序dgesv: 系统正好是奇异的
是不是迭代函数存在问题啊，谢谢了

function Newtons()
当然有问题
因为y是initial values
而你又把 y-solve(obj$J, obj$f)给了y
程序开始应该是
y1 <- y;
obj <- funs(y1);#初始值带入funs
ynew<- y1-solve(obj$J, obj$f);

你的错误来自于solve(obj$J, obj$f)
错误于drop(.Call("La_dgesv", a, as.matrix(b), tol, PACKAGE = "base")) :
  Lapack例行程序dgesv: 系统正好是奇异的

表示要迭代部分函数有误,or improper initial values.

把所求方程的一阶导数求出来（关于变量），若变量多于一个，则所求的一阶导数是一个向量，按照newton-raphson公式，最后是一个线性方程组，利用迭代法就可算出！

飘洒 发表于 2010-6-22 20:07
把所求方程的一阶导数求出来（关于变量），若变量多于一个，则所求的一阶导数是一个向量，按照newton-raphson公式，最后是一个线性方程组，利用迭代法就可算出！

谢谢回复，一阶偏倒的J和f已经算出来，现在的问题就是牛顿迭代函数有问题

1. # [Jones2009]Introduction to Scientific Programming and Simulation Using R,1420068725
   
2. newtonraphson <- function(ftn, x0, tol = 1e-9, max.iter = 100) {
   
3.         # Newton_Raphson algorithm for solving ftn(x)[1] == 0
   
4.         # we assume that ftn is a function of a single variable that returns
   
5.         # the function value and the first derivative as a vector of length 2
   
6.         #
   
7.         # x0 is the initial guess at the root
   
8.         # the algorithm terminates when the function value is within distance
   
9.         # tol of 0, or the number of iterations exceeds max.iter
   
10.         # initialise
    
11.         x <- x0
    
12.         fx <- ftn(x)
    
13.         iter <-  0
    
14.         # continue iterating until stopping conditions are met
    
15.         while ((abs(fx[1]) > tol) && (iter < max.iter)) {
    
16.                 x <- x - fx[1]/fx[2]
    
17.                 fx <- ftn(x)
    
18.                 iter <-  iter + 1
    
19.                 cat("At iteration", iter, "value of x is:", x, "\n")
    
20.         }
    
21. 
    
22.         # output depends on success of algorithm
    
23.         if (abs(fx[1]) > tol) {
    
24.                 cat("Algorithm failed to converge\n")
    
25.                 return(NULL)
    
26.         } else {
    
27.                 cat("Algorithm converged\n")
    
28.                 return(x)
    
29.         }
    
30. }
    
31. 
    
32. ftn4 <- function(x) {
    
33.     fx <- log(x) - exp(-x)
    
34.     dfx <- 1/x + exp(-x)
    
35.     return(c(fx, dfx))
    
36. }
    
37. #newtonraphson(ftn4, 2, 1e-06)

复制代码

5# epoh

谢谢，不过牛顿迭代的源程序这样写的，是步步逼近的原则。

通常把F和Jacobian Matrix分别定义,
方便在多个联立方程时除错
底下是简单两个联立方程之求解:
#function values
fun = function (x)
{
f = array(0,c(2,1))
f[1] = x[1] + 2*x[2] - 2
f[2] = x[1]^2 + 4*x[2]^2 - 4
f
}
############
#derivative values,Jacobian Matrix 2x2
Jac = function (x)
{
j = rbind(c(1,2),c(2*x[1],8*x[2]))
j
}
#############
NR = function (f, J, x0, tol, max.it)
{
xold = x0
xnew = xold - solve(Jac(xold)) %*% fun(xold)

for (k in 1:max.it)
{
xold = xnew
niter = k
xnew = xold - (solve(Jac(xold)) %*% fun(xold))
if((svd(fun(xnew))$d < tol) |  max(abs(xold-xnew))/max(abs(xnew))<tol | (niter==max.it) ) break   
xold=xnew
}                                                                                                                       
list(r=xnew,niter=niter)
}  
############
result=NR (fun, Jac, x0=c(1,2), tol=1.0e-06, max.it=20)
result
$r
             [,1]
[1,] -6.27144e-09
[2,]  1.00000e+00

$niter
[1] 4

epoh 发表于 2010-6-24 10:12
function Newtons()
当然有问题
因为y是initial values
而你又把 y-solve(obj$J, obj$f)给了y
程序开始应该是
y1

这个回复并不是最佳答案，但是您关于初值以及迭代函数的问题，从普适性来说是正确的，感谢您的热心回复。

依code,我解决了六个方程的管网分析,谢谢.

请问，能告诉一下牛顿法迭代求零点的程序吗？谢谢