Newton Raphson method in Matlab?
The Newtons-Raphsons method is easy to use in Mathematica, but it seems a little tricky in Matlab. I don't understand if I can pass a function to a function and how to use a 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 looks like you cant get neither functions nor functions defined in the file, but I could be wrong.
source to share
You can use the 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 the function descriptors of both the function and its derivative. This derivative can be acquired in several ways: manual differentiation, symbolic differentiation, or automatic differentiation. You can also do differentiation numerically, but this is slow and requires a modified implementation. So I will assume that you have computed 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);
source to share
It is impossible to algebraically take derivatives from functions or functions defined in m files. You would have to do this numerically by evaluating the function at multiple points and approximating the derivative.
What you probably want to do is differentiate symbolic equations and you need Symbolic Math Toolbox to do this. Here's an example of finding the 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 is xRoot
approaching the true root value (which is 2) after just a few iterations. You can also place the function evaluation in a while loop with a condition that checks how big the difference between each new waste and the previous guess is, stopping when the difference is small enough (i.e. the root is 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
source to share
% Friday June 07 by Ehsan Behnam.
% b) Newton 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) = [];
source to share