Newtons-Raphsons method is easy to implement in Mathematica but in Matlab it seems a bit difficult. I don't get if I can pass a function to a function and how to use the derivative as a function.
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]]]
It doesn't seem like you can derive neither function-handles nor functions defined in a file but I might be wrong.
You could use an implementation like this:
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
And pass it function handles of both the function and its derivative. This derivative is possible to acquire by some different methods: manual differentiation, symbolic differentiation or automatic differentiation. You can also perform the differentiation numerically, but this is both slow and requires you to use a modified implementation. So I will assume you have calculated the derivative in any suitable way. Then you can call the code:
f = @(x)((x-4).^2-4);
dfdx = @(x)(2.*(x-4));
x0 = 1;
xRoot = newton(@f,@dfdx,x0,1e-10);
There's no way to algebraically take derivatives of function handles or functions defined in m-files. You would have to do this numerically by evaluating the function at a number of points and approximating the derivative.
What you're probably wanting to do is differentiation of symbolic equations, and you need the Symbolic Math Toolbox for that. Here's an example of finding a root using the 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
You can see that the value of xRoot comes close to the value of the true root (which is 2) after just a couple of iterations. You could also place the function evaluation in a while loop with a condition that checks how big a difference there is between each new guess and the previous guess, stopping when that difference is sufficiently small (i.e. the root has been found):
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
% 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) = [];
