求助 关于Newton-Raphson

xilanhua

求助，有没有高人知道怎么在Matlab里面应用Newton-Raphone关于求解非线性方程组的呀。
Help里很模糊。
多谢了

mdeng1985

基本思想：
Newton-Raphson方法的优点是简单快速，缺点是需要目标函数的一阶导数df/dx,因此这种方法适用于目标函数导数可求的情况下。
以单变量的函数为例，将f在xi点taylor展开：
f(xi+1)=f(xi) + f’(xi) ( xi+1- xi ) +O(xi+1- xi)^2
略去高阶项，如果xi+1是方程的根的话，则
xi+1= xi- f(xi)/ f’(xi)
这样就可以由上一个迭代点推算出下一个迭代点,直到abs(xi+1- xi)满足一定精度为止
以上就是Newton-Raphson的基本思想

算法实现：
1，        给定一个初始估计x0
2，        计算delta x=- f(xi)/ f’(xi)
3，        xi+1= xi+detla x，
4，        判断abs(xi+1- xi)是否小于epilson(给定的精度)，如果满足则xi+1就是求得的非线性方程的根，如果不满足则继续2的工作，直到满足为止

多变量的情况和单变量的情况原理一样，只用把第2步换成计算多变量的偏导数，也就是jacobian矩阵
matlab实现
方法1，根据自己的理解，但不一定对
options=optimset('NonlEqnAlgorithm','gn ')
options=optimset(options,'Jacobian','on')
定义目标函数
function [f,g]=myfun(x)
f=目标函数的解析式子
g=梯度函数的解析式子

方法2，保守做法，把下面的代码保存为m文件，调用即可
function root = newtonRaphson2(func,x,tol)
% Newton-Raphson method of finding a root of simultaneous
% equations fi(x1,x2,...,xn) = 0, i = 1,2,...,n.
% USAGE: root = newtonRaphson2(func,x,tol)
% INPUT:
% func = handle of function that returns[f1,f2,...,fn].
% x = starting solution vector [x1,x2,...,xn].
% tol = error tolerance (default is 1.0e4*eps).
% OUTPUT:
% root = solution vector.
if nargin == 2; tol = 1.0e4*eps; end
if size(x,1) == 1; x = x'; end % x must be column vector
for i = 1:30
[jac,f0] = jacobian(func,x);
if sqrt(dot(f0,f0)/length(x)) < tol
root = x; return
end
dx = jac\(-f0);
x = x + dx;
if sqrt(dot(dx,dx)/length(x)) < tol*max(abs(x),1.0)
root = x; return
end
end
error('Too many iterations')
function [jac,f0] = jacobian(func,x)
% Returns the Jacobian matrix and f(x).
h = 1.0e-4;
n = length(x);
jac = zeros(n);
f0 = feval(func,x);
for i =1:n
temp = x(i);
x(i) = temp + h;
f1 = feval(func,x);
x(i) = temp;
jac(:,i) = (f1 - f0)/h;
end

[ 本帖最后由 mdeng1985 于 2007-10-13 01:34 编辑 ]

xilanhua

太感谢mdeng1985这位朋友了，虽然还没有具体实现，但是对我肯定是帮助的。太感谢了。

淘汰郎

请教三楼，方法2运行时候出现这样一个错误是怎么回事？
Input argument "x" is undefined.

Error in ==> newtonRaphson2 at 12
if size(x,1) == 1; x = x; end % x must be column vecto

lin2009

4# 淘汰郎
提示信息很简单啊。
需要一个列向量，可能你的x为其它的数据结构。如行向量或矩阵等。

另外，需要把代码和出错信息都粘上来才意义。