Newtons-Raphsons 法は Mathematica では簡単に実装できますが、Matlab では少し難しいようです。関数を関数に渡すことができるかどうか、および導関数を関数として使用する方法がわかりません。
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]]]
ファイルで定義された関数ハンドルも関数も派生できないようですが、私は間違っているかもしれません。
次のような実装を使用できます。
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);
関数ハンドルまたは m ファイルで定義された関数の導関数を代数的に取得する方法はありません。いくつかの点で関数を評価し、導関数を近似することにより、これを数値的に行う必要があります。
おそらくやりたいことは、シンボリック方程式の微分であり、そのためにはSymbolic Math Toolboxが必要です。ニュートン・ラフソン法を使用して根を見つける例を次に示します。
>> 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
% 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) = [];