這篇文章給出不動點及牛頓法迭代是完整實現步驟，沒有考慮運行效率及內存損耗，僅供參考。
一、舉例計算
求下列方程的實根
(1) . x^2 - 3x + 2 - exp(x) = 0
(2). x^3 + 2x^2 +10x - 20 = 0
設計一種不動點及牛頓法迭代法，使迭代序列收斂，且計算到x(k) - x(k-1)< 10-8為止，結果輸出小數點 12 位及以上，通過迭代次數及精度等來比較方法的優劣。
二、不動點及牛頓法迭代實現
2.1 、實現環境
MATLAB 2018b
2.2 、實現原理
不動點迭代法的實現:通過合理構造j(x) 函數，滿足不動點迭代局部收斂的條件：即在真值 x* 的鄰域附近使得|j'(x)| < 1，在程序中首先需要大概計算下構造函數的導數，看是否滿足條件，以其為后面迭代函數，而在不動點迭代函數MATLAB實現中，編寫了一個函數myALG(x0, func, x_real, tol, MaxIter) ，輸入參數有初值、構造的j(x) 函數、經過vpasolve 計算得出的真實值、誤差容忍值及最大迭代次數參數。將這些參數傳輸至函數中后，函數通過定義中間變量來存儲迭代值，將初值送入至構造的迭代函數后的輸出賦值給中間變量，然后依據設置的誤差容忍程度來進行循環。
而牛頓法迭代的實現需要重新構造j(x) 函數：j(x) = x - f(x)/f'(x) ，需要對所要計算的方程進行求導等進行構造，而在 MATLAB的實現也是和不動點迭代法類似，都是函數通過定義中間變量來存儲迭代值，將初值送入至構造的迭代函數后的輸出賦值給中間變量，然后依據設置的誤差容忍程度來進行循環。兩種方法的具體實現見附錄代碼。
2.3 、數據測試
方程 1 的實現： x^2 - 3x + 2 - exp(x) = 0
不動點迭代構造j(x) 函數：j(x) = x^2 / 3 + 2 / 3 - exp(x) / 3
牛頓法迭代構造j(x) 函數：j(x) = x - (x^2 - 3x + 2 - exp(x) ) /(2x - 3 - exp(x))初始化myALG(2.0,@ myFunc, real_value,1e-8,20) ，之后通過 MATLAB 進行測試，經過vpasolve函數計算方程1，得到方程1解為x* =0.25753028543986076045536730493724，之后進行不動點及牛頓法迭代測試。
不動點迭代得到如下結果:
牛頓法迭代得到如下結果:
方程 2 的實現： x3 + 2x2 +10x - 20 = 0
不動點迭代構造j(x) 函數：j(x) = x + (x^3 + 2x^2 +10x - 20) × (-1/ 20)
牛頓法迭代構造j(x) 函數：j(x) = x - (x^3 + 2x^2 +10x - 20) /(3x^2 + 4x +10)
初始化myALG(1.3@ myFunc, real _ value,1e - 8,15) ，之后通過 MATLAB 進行測試， 經過 vpasolve 函數計算方程 2，得到方程 2有幾個不同解， 以解x* =1.3688081078213726352274143300213 為例來進行測試。
不動點迭代得到如下結果:
牛頓法迭代得到如下結果:
2.4 、數據分析
通過構造了合適的不動點及牛頓法迭代，實現了方程 1 與方程 2 的求解，並使之與 xk - xk -1 < 10-8 及結果輸出小數點 12 位及以上等要求。經過測試可以看出在一定條件下牛頓法收斂的速度要比不動點迭代法的收斂速度要快，甚至相差一定數量級，這與數值分析課本中對牛頓法及不動點迭代法收斂速度相對應。
附錄
不動點迭代代碼：
math_test_BDD.m
1 % syms x 2 % real_value=vpasolve(x^2-3*x+2-exp(x) == 0, x) %不動點迭代方程式1 3 % cal_value=myALG(2.0,@myFunc,real_value,1e-8,15); 4 % 5 % if abs(cal_value - real_value)< 1e-8 6 % fprintf('abs(cal_value - real_value)*10^8=%f\n',abs(cal_value - real_value)*1e8); 7 % end 8 % abs(cal_value - real_value) 9 10 11 syms x 12 real_value=vpasolve(x^3 + 2*x^2 + 10*x - 20 == 0, x) %不動點迭代方程式2 13 real_value_r = real_value(1) 14 cal_value=myALG_BDD(1.3,@myFunc_BDD,real_value_r,1e-8,15); 15 16 if abs(cal_value - real_value_r)< 1e-8 17 fprintf('abs(cal_value - real_value_r)*10^8=%f\n',abs(cal_value - real_value_r)*1e8); 18 end 19 abs(cal_value - real_value_r)
myALG_BDD.m
1 function [y] = myALG_BDD(x0,func,x_real,tol,MaxIter) 2 3 xn = x0; 4 fprintf('Iter 0: %16.14f\n',x0); 5 xnp1 = func(xn); 6 fprintf('Iter 1: %16.14f\n',xnp1); 7 %criterion = abs(xnp1-xn); 8 criterion = abs(xnp1-x_real); 9 xn = xnp1; 10 Iter = 1; 11 %while(criterion>tol) 12 while(abs(criterion)>tol) 13 xnp1 = func(xn); 14 %criterion = abs(xnp1-xn); 15 criterion = abs(xnp1-x_real); 16 xn = xnp1; 17 Iter = Iter + 1; 18 fprintf('Iter %2.0d: %16.14f\n',Iter,xnp1); 19 if Iter>=MaxIter 20 break; 21 end 22 end 23 y = xnp1; 24 end
myFunc_BDD.m
1 function [y] = myFunc_BDD(x) 2 3 %y = x^2-3*x+2-exp(x); 4 %y = x^2/3+2/3-exp(x)/3; %不動點迭代方程式1 5 6 %y = x^3 + 2*x^2 + 10*x - 20 y = -x^3/10 - x^2/5 + 2 7 y = x + (x^3 + 2*x^2 + 10*x - 20)*(-1/20); %不動點迭代方程式2 8 end
牛頓法迭代代碼：
math_test_ND.m
1 % syms x 2 % real_value=vpasolve(x^2-3*x+2-exp(x) == 0, x) %不動點迭代方程式1 3 % cal_value=myALG(2.0,@myFunc,real_value,1e-8,15); 4 % 5 % if abs(cal_value - real_value)< 1e-8 6 % fprintf('abs(cal_value - real_value)*10^8=%f\n',abs(cal_value - real_value)*1e8); 7 % end 8 % abs(cal_value - real_value) 9 10 11 % syms x 12 % real_value=vpasolve(x^3 + 2*x^2 + 10*x - 20 == 0, x) %不動點迭代方程式2 13 % real_value_r = real_value(1) 14 % cal_value=myALG(1.3,@myFunc,real_value_r,1e-8,15); 15 % 16 % if abs(cal_value - real_value_r)< 1e-8 17 % fprintf('abs(cal_value - real_value_r)*10^8=%f\n',abs(cal_value - real_value_r)*1e8); 18 % end 19 % abs(cal_value - real_value_r) 20 21 22 23 % syms x 24 % real_value=vpasolve(x^2-3*x+2-exp(x) == 0, x) %牛頓法迭代1 25 % real_value_r = real_value 26 % cal_value=myALG(2,@myFunc,real_value_r,1e-8,15); 27 % 28 % if abs(cal_value - real_value_r)< 1e-8 29 % fprintf('abs(cal_value - real_value_r)*10^8=%f\n',abs(cal_value - real_value_r)*1e8); 30 % end 31 % abs(cal_value - real_value_r) 32 33 34 syms x 35 real_value=vpasolve(x^3 + 2*x^2 + 10*x - 20 == 0, x) %牛頓法迭代2 36 real_value_r = real_value(1) 37 cal_value=myALG_ND(1.3,@myFunc_ND,real_value_r,1e-8,15); 38 39 if abs(cal_value - real_value_r)< 1e-8 40 fprintf('abs(cal_value - real_value_r)*10^8=%f\n',abs(cal_value - real_value_r)*1e8); 41 end 42 abs(cal_value - real_value_r)
myALG_ND.m
1 function [y] = myALG_ND(x0,func,x_real,tol,MaxIter) 2 3 xn = x0; 4 fprintf('Iter 0: %16.14f\n',x0); 5 xnp1 = func(xn); 6 fprintf('Iter 1: %16.14f\n',xnp1); 7 %criterion = abs(xnp1-xn); 8 criterion = abs(xnp1-x_real); 9 xn = xnp1; 10 Iter = 1; 11 %while(criterion>tol) 12 while(abs(criterion)>tol) 13 xnp1 = func(xn); 14 %criterion = abs(xnp1-xn); 15 criterion = abs(xnp1-x_real); 16 xn = xnp1; 17 Iter = Iter + 1; 18 fprintf('Iter %2.0d: %16.14f\n',Iter,xnp1); 19 if Iter>=MaxIter 20 break; 21 end 22 end 23 y = xnp1; 24 end
myFunc_ND.m
1 function [y] = myFunc_ND(x) 2 3 4 5 %y = x^2-3*x+2-exp(x); 6 %y = x^2/3+2/3-exp(x)/3; %不動點迭代方程式1 7 8 %y = x^3 + 2*x^2 + 10*x - 20 y = -x^3/10 - x^2/5 + 2 9 %y = x + (x^3 + 2*x^2 + 10*x - 20)*(-1/20); %不動點迭代方程2 10 11 %y = x-(x^2-3*x+2-exp(x))/(2*x-3-exp(x)); %牛頓法迭代1 12 13 y = x - (x^3 + 2*x^2 + 10*x - 20)/(3*x^2+4*x+10); %牛頓法迭代2 14 15 16 end