% Newton法サンプル
clear all;
clc;
close all;

% 諸設定
x = 0.1;        % これを色々変えてみると楽しい
                % お薦め値： -2.2（快適) 1.6(やや快適) 1.7(迷走) 0.1(迷走)
tol = 1e-6;     % |f(x)|<tol となったら停止（停止則）
format long;    % MATLABの表示桁数を倍精度相当にするおまじない

% Newton 法の反復
y = funcy(x);   % y(x) を2回計算しないためのちまちました工夫
xl = x;         % xl = {近似解のリスト}
while 1     % 無限ループ
    x = x -  y / derivy(x);     % Newton 反復式
    x                           % x を表示
    xl = [xl x];                % {近似解のリスト} を更新
    y = funcy(x);               % 新しい x における残差 f(x)
    if abs(y) < tol
        break;                  % 残差が許容誤差以下になったら脱出
    end
end

% サービス：関数を図示
xx = linspace(-2, 3, 100);      % 図示する範囲の分点を指定
yy = funcy(xx);     % xx に対して一気に y(xx) を計算（MATLABの機能）
plot(xx, yy);       % y(x) を描画
hold on;            % 図を消さないで上書きして，というお願い
plot(xx, zeros(length(xx)));    % 実軸を描画
plot(xl, 0, 'go');     % {近似解のリスト} を描画

% 関数の定義：x^3 - 2x^2 - x + 4
function y = funcy(x)
    y = x.^3 - 2.0 * x.^2 - x + 4.0;    % x.^3 というのはベクトル x の各要素を 3 乗する，という意味
                                        % こう書いておくとこの関数はスカラーにもベクトルにも使える
end

% 関数の微分の定義（人間が手動で・・・）：3x^2 - 4x - 1
%   一般には Jacobian を人間が計算してプログラミングするので大変
%   微分を「差分」で近似する「数値 Jacobian」という概念もある
%   微分値を，計算グラフに沿ってアルゴリズミックに計算する「高速自動微分」という概念もある
%   （近年深層学習で使われている後退誤差伝播法(back propagation)は実は自動微分）
function dy = derivy(x)
    dy = 3.0 * x^2 - 4.0 * x - 1.0;
end