2011-04-12 38 views
1

Newtons-Raphsons方法在Mathematica中很容易实现,但在Matlab中似乎有点困难。我不明白我是否可以将函数传递给函数,以及如何将函数用作函数。Matlab中的Newton Raphsons方法?

newtonRaphson[f_, n_, guess_] := 
If[n == 0, guess, newtonRaphson[f, n - 1, guess - f[guess]/f'[guess]]] 
newtonRaphsonOptimize[f_, n_, guess_] := 
If[n == 0, guess, 
    newtonRaphsonOptimize[f, n - 1, guess - f'[guess]/f''[guess]]] 

它似乎并不像你可以派生函数句柄或函数定义的文件,但我可能是错的。

回答

2

代数方式没有办法在m文件中定义的函数句柄或函数的衍生物。您必须通过在多个点处评估函数并近似导数来以的数值计算

你可能想要做的是differentiation of symbolic equations,你需要这个Symbolic Math Toolbox。下面是一个使用Newton-Raphson method找到根的例子:

>> syms x   %# Create a symbolic variable x 
>> f = (x-4)^2-4; %# Create a function of x to find a root of 
>> xRoot = 1;  %# Initial guess for the root 
>> g = x-f/diff(f); %# Create a Newton-Raphson approximation function 
>> xRoot = subs(g,'x',xRoot) %# Evaluate the function at the initial guess 

xRoot = 

    1.8333 

>> xRoot = subs(g,'x',xRoot) %# Evaluate the function at the refined guess 

xRoot = 

    1.9936 

>> xRoot = subs(g,'x',xRoot) %# Evaluate the function at the refined guess 

xRoot = 

    2.0000 

你可以看到的xRoot值接近真正的根本价值后只是一对夫妇的迭代(其中2个是)。您还可以将函数评估放置在while循环中,条件是检查每个新猜测与之前的猜测之间有多大的差异,当差异足够小时即停止(即已找到根目录)时停止:

xRoot = 1;      %# Initial guess 
xNew = subs(g,'x',xRoot);  %# Refined guess 
while abs(xNew-xRoot) > 1e-10 %# Loop while they differ by more than 1e-10 
    xRoot = xNew;    %# Update the old guess 
    xNew = subs(g,'x',xRoot); %# Update the new guess 
end 
xRoot = xNew;     %# Update the final value for the root 
+0

感谢,但我不明白如何使用,要么。我需要以某种方式将函数传递给newtons raphsons。以及如何评估做差异(f)后? – jon 2011-04-12 18:07:56

+0

而不是将整个函数传递给另一个函数,为什么不考虑传递结果? – 2013-06-08 07:38:05

6

你可以使用这样的实现:

function x = newton(f,dfdx,x0,tolerance) 
err = Inf; 
x = x0; 
while abs(err) > tolerance 
    xPrev = x; 
    x = xPrev - f(xPrev)./dfdx(xPrev); 
    % stop criterion: (f(x) - 0) < tolerance 
    err = f(x); % 
    % stop criterion: change of x < tolerance 
    % err = x - xPrev; 
end 

而且它传递函数及其导数两者的函数处理。这种衍生物有可能通过一些不同的方法获得:手动分化,符号分化或自动分化。您也可以用数字来执行区分,但这样做速度很慢,并且需要您使用修改的实现。所以我会假设你已经以任何适当的方式计算了衍生物。然后,你可以调用代码:

f = @(x)((x-4).^2-4); 
dfdx = @(x)(2.*(x-4)); 
x0 = 1; 
xRoot = newton(@f,@dfdx,x0,1e-10); 
0
% Friday June 07 by Ehsan Behnam. 
% b) Newton's method implemented in MATLAB. 
% INPUT:1) "fx" is the equation string of the interest. The user 
% may input any string but it should be constructable as a "sym" object. 
% 2) x0 is the initial point. 
% 3) intrvl is the interval of interest to find the roots. 
% returns "rt" a vector containing all of the roots for eq = 0 
% on the given interval and also the number of iterations to 
% find these roots. This may be useful to find out the convergence rate 
% and to compare with other methods (e.g. Bisection method). 
% 
function [rt iter_arr] = newton_raphson(fx, x, intrvl) 
n_seeds = 10; %number of initial guesses! 
x0 = linspace(intrvl(1), intrvl(2), n_seeds); 
rt = zeros(1, n_seeds); 

% An array that keeps the number of required iterations. 
iter_arr = zeros(1, n_seeds); 
n_rt = 0; 

% Since sometimes we may not converge "max_iter" is set. 
max_iter = 100; 

% A threshold for distinguishing roots coming from different seeds. 
thresh = 0.001; 

for i = 1:length(x0) 
    iter = 0; 
    eq = sym(fx); 
    max_error = 10^(-12); 
    df = diff(eq); 
    err = Inf; 
    x_this = x0(i); 
    while (abs(err) > max_error) 
     iter = iter + 1; 
     x_prev = x_this; 

     % Iterative process for solving the equation. 
     x_this = x_prev - subs(fx, x, x_prev)/subs(df, x, x_prev); 
     err = subs(fx, x, x_this); 
     if (iter >= max_iter) 
      break; 
     end 
    end 
    if (abs(err) < max_error) 
     % Many guesses will result in the same root. 
     % So we check if the found root is new 
     isNew = true; 
     if (x_this >= intrvl(1) && x_this <= intrvl(2)) 
      for j = 1:n_rt 
       if (abs(x_this - rt(j)) < thresh) 
        isNew = false; 
        break; 
       end 
      end 
      if (isNew) 
       n_rt = n_rt + 1; 
       rt(n_rt) = x_this; 
       iter_arr(n_rt) = iter; 
      end 
     end 
    end   
end 
rt(n_rt + 1:end) = []; 
iter_arr(n_rt + 1:end) = [];