一、非線性方程求根
通過以下問題學習此知識點：
現在你想買一套300萬元的房子，首付40%，貸款20年，等額本息，已知月還款額為1.2萬元，求貸款月利率為多少？
(1) 編寫結合牛頓下山法和割線法的綜合疊代方法求解函數，調用後求解；
(2) 使用steffenson法求解。
1、牛頓疊代法
- 又稱為牛頓-拉弗森方法（Newton-Raphson method），單變量下又稱為切線法。它是一種在實數域和複數域上近似求解方程的方法。方法使用函數f (x)的泰勒級數的前面幾項來尋找方程f (x) = 0的根。用牛頓疊代法解非線性方程，是把非線性方程f(x) = 0線性化的一種近似方法。
-
牛頓疊代法的本質是一種線性化方法，其基本思想是将非線性方程f(x)=0逐漸歸結為某種線性方程來求解。
設方程f(x)=0有近似根，将函數f(x)在點處展開，則有
于是，方程f(x)=0就可以近似表示為
非線性方程組求解-MatLab
對該方程求解，得到
x k + 1 = x k − f ( x k ) f ′ ( x k ) ( k = 0 , 1 , 2 , 3... ) x_{k+1}=x_k-\frac{f(x_k)}{f^{'}(x_k)} \quad\quad (k=0,1,2,3...) xk+1=xk−f′(xk)f(xk)(k=0,1,2,3...)
該疊代公式即為 牛頓疊代公式。
x k + 1 x_{k+1} xk+1是 y = f ( x ) y=f(x) y=f(x)的切線在點 ( x k , f ( x k ) ) (x_k,f(x_k)) (xk,f(xk))處與 x x x軸交點的橫坐标。
算法流程圖：
對牛頓疊代法更加直覺的了解
2、牛頓下山法
在牛頓疊代法中，有時候會出現如下疊代回退的情況，導緻出現死循環
當選取初值有困難時,可改用如下疊代格式,以擴大初值的選取範圍,
x n + 1 = x n − λ f ( x k ) f ′ ( x k ) x_{n+1}=x_n-\lambda\frac{f(x_k)}{f^{'}(x_k)} xn+1=xn−λf′(xk)f(xk)
其中λ稱為下山因子,λ選取應滿足單調性條件
∣ f ( x n + 1 ∣ < ∣ f ( x n ) ∣ |f(x_{n+1}|<|f(x_n)| ∣f(xn+1∣<∣f(xn)∣
這樣的算法稱下山法.将下山法和牛頓法結合起來使用的方法,稱為牛頓下山法.
下山因子λ的選擇是逐漸探索的過程.從λ=1開始反複将λ減半進行試算,如果能定出λ使單調性條件成立則“下山成功”.與此相反,如果找不到使單調性條件成立的λ,則“下山失敗”.此時需另選初值x0重算
3、割線法引入
割線法又稱為弦截法法
-
為什麼要引入割線法？
在牛頓下山法中我們需要計算被切割點的倒數，在多次疊代情況下，調用MatLab自帶的求導函數會占用很多計算資源，是以，我們采用向後有限差分來近似導數計算，什麼意思？
上圖：
非線性方程組求解-MatLab 在割線法中，省略了求導步驟，用 x k x_k xk和 x k − 1 x_{k-1} xk−1兩點連成的直線與x軸的交點近似代替了牛頓法中 x k + 1 x_{k+1} xk+1的求解
由此，減少了計算資源的消耗
-
需要注意
弦截法和牛頓疊代法都是線性化方法,牛頓疊代法在計算 x n + 1 x_{n+1} xn+1時隻用到前一步的值 x n x_n xn,而弦截法用到前面兩步的結果 x n x_n xn和 x n + 1 x_{n+1} xn+1,是以使用割線法必須先給出兩個初值值 x 0 , x 1 x_0,x_{1} x0,x1.
牛頓下山法和割線法結合算法如下：
function [x,iter,X]=newton_secant(fun,x0,x1,eps,maxiter)
% 牛頓下山+割線法求解非線性方程的根
% 輸入參數：
% ---fun：疊代函數
% ---x0,x1：初始疊代點
% ---eps：精度要求，預設值為1e-6
% ---maxiter：最大疊代次數，預設值為1e4
% 輸出參數：
% ---x：非線性方程的近似根
% ---iter：疊代次數
% ---X：每一步疊代的結果
if nargin<3,error('輸入參數至少需要3個！'),end
if nargin<4|isempty(eps),eps=1e-6;end
if nargin<5|isempty(maxiter),maxiter=1e4;end
f0 = feval(fun,x0); % 計算x0處的函數值
f1 = feval(fun,x1); % 計算x1處的函數值
k=0;err=1;%k統計疊代次數，err表示疊代精度
while abs(err)>eps
lambda=1;%下山因子
x2=x1-lambda*f1*(x1-x0)/(f1-f0);
f2 = feval(fun,x2);
while abs(f2)>=abs(f1)
lambda=lambda/2; % 更新lambda
x2=x1-lambda*f1*(x1-x0)/(f1-f0); % 牛頓下山疊代,割線法代替求導
f2=feval(fun,x2);
end
x0=x1;x1=x2;
f0=f1;f1=f2;
err=x1-x0;%更新疊代精度
k=k+1;
X(k)=x2;
end
if k>=maxiter
error('疊代次數超限，疊代失敗！')
end
x=x2;iter=k;X=X';
end
解決第一個問題
- 每月還款數額計算公式如圖：
非線性方程組求解-MatLab P:貸款本金
R:月利率
N:還款期數
月利率 = 年利率/12
- 是以我們有
180 × x ( 1 + x ) 240 ( 1 + x ) 240 − 1 − 1.2 = 0 180\times\frac{x(1+x)^{240}}{(1+x)^{240}-1}-1.2=0 180×(1+x)240−1x(1+x)240−1.2=0
fun=@(x)180*(x*(1+x)^240)/((1+x)^240-1)-1.2;
指令行調用函數解方程：
[x,iter,X]=newton_secant(fun,0.007,0.006)
輸出：
x =
0.00426762528564472
iter =
4
X =
0.00437256674949125
0.0042721389151354
0.00426763795820826
0.00426762528564472
可以看到：
月利率約為 0.427 % 0.427\% 0.427%
4、steffenson法
Steffeson法的導數近似方法：
f ′ ( x k ) ≈ f ( x k + f ( x k ) ) − f ( x k ) f ( x k ) f^{'}(x_k)\approx\frac{f(x_k+f(x_k))-f(x_k)}{f(x_k)} f′(xk)≈f(xk)f(xk+f(xk))−f(xk)
則疊代公式為：
x k + 1 = x k − f 2 ( x k ) f ( x k + f ( x k ) ) − f ( x k ) ( k = 1 , 2 , 3... ) x_{k+1}=x_k-\frac{f^{2}(x_k)}{f(x_k+f(x_k))-f(x_k)}\quad\quad (k=1,2,3...) xk+1=xk−f(xk+f(xk))−f(xk)f2(xk)(k=1,2,3...)
廢話少說，上代碼：
function [x,iter,X] = steffenson(fun,x0,eps,maxiter)
% Steffenson法求解非線性方程的根
% 輸入參數：
% ---fun：疊代函數
% ---x0：初始疊代點
% ---eps：精度要求，預設值為1e-6
% ---maxiter：最大疊代次數，預設值為1e4
% 輸出參數：
% ---x：非線性方程的近似根
% ---iter：疊代次數
% ---X：每一步疊代的結果
if nargin<2,error('輸入參數至少需要2個！'),end
if nargin<3|isempty(eps),eps=1e-6;end
if nargin<4|isempty(maxiter),maxiter=1e4;end
k=0;err=1;
while abs(err)>eps;
k=k+1;
f0 = feval(fun,x0); % 計算x0處的函數值
x1=x0-f0^2/(feval(fun,x0+f0)-f0); % Steffenson法疊代公式
err=x1-x0;
x0 = x1; % 更新x0數值
X(k)=x1;
end
x=x1;iter=k;X=X';
指令行調用函數解方程：
[x,iter,X]=steffenson(fun,0.007,0.006)
輸出：
x =
0.00426762553393485
iter =
6
X =
0.005042607730931
0.0045032072675304
0.00433124100792348
0.00427709276934277
0.00426790263879553
0.00426762553393485
可以看到：
月利率約為 0.427 % 0.427\% 0.427%，但是疊代次數為6次
二、非線性方程組求解
通過以下問題學習此知識點：
用牛頓法求解二進制方程組的根
x 2 c o s 2 x + y 2 s i n ( 2 y ) = 1 x 3 + y 3 − 6 c o s ( 2 x y ) = − 1 \begin{array}{cc} x^2 \mathrm{cos2x}+y^2 \mathrm{sin(2y)}=1 & \\ x^3 +y^3 -6\mathrm{cos(2xy)}=-1 & \end{array} x2cos2x+y2sin(2y)=1x3+y3−6cos(2xy)=−1
1、牛頓法解方程組
非線性方程組的求法有很多，此處僅對使用牛頓法求解非線性方程組的根進行學習。
将多元向量函數F(x) 在點處展開
F ( x ) ≈ F ( x ( k ) ) + F ′ ( x ( k ) ) ( x − x ( k ) ) ) F\left(x\right)\approx F\left(x^{\left(k\right)} \right)+F^{\prime } \left(\left.x^{\left(k\right)} \right)\left(x-x^{\left(k\right)} \right)\right) F(x)≈F(x(k))+F′(x(k))(x−x(k)))
其中，是F(x)的Jacobi矩陣
是以，可以得到求解非線性方程組的疊代方程
x ( k + 1 ) = x ( k ) − [ F ′ ( x ( k ) ) ] − 1 F ( x ( k ) ) x^{\left(k+1\right)} =x^{\left(k\right)} -{\left\lbrack F^{\prime } \left(x^{\left(k\right)} \right)\right\rbrack }^{-1} F\left(x^{\left(k\right)} \right) x(k+1)=x(k)−[F′(x(k))]−1F(x(k))
---------流程圖--------
上代碼：
function [x,iter,X]=newtong(fun,x0,eps,maxiter)
% Newton法求解非線性方程組的根
% 輸入參數：
% ---fun：疊代函數
% ---x0：初始疊代點向量
% ---eps：精度要求，預設值為1e-6
% ---maxiter：最大疊代次數，預設值為1e4
% 輸出參數：
% ---x：非線性方程的近似解向量
% ---iter：疊代次數
% ---X：每一步疊代的結果
if nargin<2,error('輸入參數至少需要2個！'),end
if nargin<3|isempty(eps),eps=1e-6;end
if nargin<4|isempty(maxiter),maxiter=1e4;end
k=0;err=1;
while err>eps
k=k+1;
[fx0,J]=feval(fun,x0); % 求函數fun的值和jacobi矩陣
x1=x0-J\fx0; % 牛頓法疊代公式
err=norm(x1-x0);
x0=x1;
X(k,:)=x1;
end
if k==maxiter
error('疊代次數超限，疊代失敗！')
end
x=x1;iter=k;
end
function [y,J]=fun(x)
% 非線性方程組
% 函數檔案描述，傳回函數值和jacobi矩陣
y=[x(1)^2*cos(2*x(1))+x(2)^2*sin(2*x(2))-1;
x(1)^3+x(2)^3-6*cos(2*x(1)*x(2))+1];
% 求Jacobi矩陣
syms xx yy; % 聲明符号變量
%J=jacobian([2*xx-yy-exp(-xx);-xx+2*yy-exp(-yy)],[xx yy]); % 求符号jacobi矩陣
J = jacobian([xx^3*cos(2*xx)+yy^2*sin(2*yy)-1,xx^3+yy^3-6*cos(2*xx*yy)+1], [xx, yy]);
xx = x(1);
yy = x(2);
J=eval(J); % 替換
end
指令行輸入：
[x,iter,X]=newtong(@fun,[5;2])
輸出：
x =
16.0945044603024
-16.1035071010038
iter =
98
X =
5.10274043149926 -0.119124207216952
5.24942127220357 2.38096639162783
5.02977093770551 -13.5938666862748
16.0396552014422 -8.16234689739867
16.3648055755991 -18.4422298393681
16.5029906700455 -16.5584665672948
……