方程求根——牛頓迭代法

這段代碼實現了牛頓切線法、簡化牛頓法和牛頓下山法這三種方程求解法，由於輸出結果較長，只以牛頓下山法為例寫一段例題

　　1.代碼

%%牛頓迭代法
%%method為-1時為牛頓切線法，method為0時為簡化牛頓法，method為1時為牛頓下山法
%%f是表達式f(x) = 0，X0是初值，epsilon是精度，interval是包含解的區間
function NM = Newton_method(f,X0,epsilon,interval,method)

Y0 = subs(f,X0);
%%作圖
t = interval(1):(interval(2)-interval(1))/1000:interval(2);
T = subs(f,t);
T1 = zeros(1,max(size(t)));
Y1 = subs(f,X0)+subs(diff(f),X0)*(t-X0);
h = figure;
set(h,'color','w');
plot(t,T,'c',t,Y1,'g',X0,Y0,'ro',t,T1,'y');
grid on;
xlabel('x shaft');ylabel('y shaft');
title('函數圖像');
hold on

x(1) = X0;
ub = 100;e = floor(abs(log(epsilon)));

if method == -1
    disp('牛頓切線法');
    for i = 2:ub
        x(i) = x(i-1)-subs(f,x(i-1))/subs(diff(f),x(i-1));
        delta = x(i)-x(i-1);
        if abs(delta) < epsilon
            break;
        end
    end
    disp('迭代次數為：');
    i-1
    disp('迭代解為：');
    NM = vpa(x,e);
    X_end = x(i);
    Y_end = subs(f,X_end);
    
    X = double([X0 X_end]);Y = double([Y0 Y_end]);
    Y2 = Y_end+subs(diff(f),X_end)*(t-X_end);
    plot(t,Y2,'b',X_end,Y_end,'mo');
    legend('T：函數圖像','Y1：初始點處切線','Y0：初始值處切點','T1：直線y=0','Y2：迭代解處的切線','Y_end：迭代解處切點');
    
    for i = 1:2
        text(X(i),Y(i),['(',num2str(X(i)),',',num2str(Y(i)),')'],'color',[0.02 0.79 0.99]);
    end
    
elseif method == 0
    disp('簡化牛頓法');
    for i = 2:ub
        x(i) = x(i-1)-subs(f,x(i-1))/subs(diff(f),x(1));
        delta = x(i)-x(i-1);
        if abs(delta) < epsilon
            break;
        end
    end
    disp('迭代次數為：');
    i-1
    disp('迭代解為：');
    NM = vpa(x,e);
    X_end = x(i);
    Y_end = subs(f,X_end);
    
    X = double([X0 X_end]);Y = double([Y0 Y_end]);
    Y2 = Y_end+subs(diff(f),X_end)*(t-X_end);
    plot(t,Y2,'b',X_end,Y_end,'mo');
    legend('T：函數圖像','Y1：初始點處切線','Y0：初始值處切點','T1：直線y=0','Y2：迭代解處的切線','Y_end：迭代解處切點');
    
    for i = 1:2
        text(X(i),Y(i),['(',num2str(X(i)),',',num2str(Y(i)),')'],'color',[0.02 0.79 0.99]);
    end
    
elseif method == 1
    disp('牛頓下山法');
    lambda = input('輸入下山因子：');
    for i = 2:ub
        x(i) = x(i-1)-lambda*subs(f,x(i-1))/subs(diff(f),x(1));
        delta = x(i)-x(i-1);
        if abs(delta) < epsilon
            break;
        end
    end
    disp('迭代次數為：');
    i-1
    disp('迭代解為：');
    NM = vpa(x,e);
    X_end = x(i);
    Y_end = subs(f,X_end);
    
    X = double([X0 X_end]);Y = double([Y0 Y_end]);
    Y2 = Y_end+subs(diff(f),X_end)*(t-X_end);
    plot(t,Y2,'b',X_end,Y_end,'mo');
    legend('T：函數圖像','Y1：初始點處切線','Y0：初始值處切點','T1：直線y=0','Y2：迭代解處的切線','Y_end：迭代解處切點');
    
    for i = 1:2
        text(X(i),Y(i),['(',num2str(X(i)),',',num2str(Y(i)),')'],'color',[0.02 0.79 0.99]);
    end
end

　　2.例子

clear all
clc
syms x;
f = x^exp(x)-1;
X0 = 0.8;
epsilon=1e-6;
interval = [0,2];
method = 1;

%%牛頓下山法
X = Newton_method(f,X0,epsilon,interval,method)

　　結果如下

牛頓下山法
輸入下山因子：0.8
迭代次數為：
ans =
    21
迭代解為：
X =
[ 0.8, 1.025142203855, 0.9839179688712, 1.008330092139, 0.995101056432, 1.002691649998, 0.9984619264687, 1.000859946529, 0.99951321039, 1.000273650001, 0.9998455622605, 1.000086966616, 0.9999509665425, 1.000027626624, 0.9999844283422, 1.000008774959, 0.9999950545033, 1.000002787045, 0.9999984292923, 1.000000885191, 0.9999995011337, 1.000000281144]

　　

由於迭代函數原因，圖象上的數據顯示出現了遮擋，這一部分的代碼以后再進行優化