MatLab非線性方程組求解--十安辰

一、非線性方程求根

通過(guò)以下問(wèn)題學(xué)習(xí)此知識(shí)點(diǎn):

現(xiàn)在你想買(mǎi)一套300萬(wàn)元的房子,首付40%,貸款20年,等額本息,已知月還款額為1.2萬(wàn)元,求貸款月利率為多少?

(1) 編寫(xiě)結(jié)合牛頓下山法和割線法的綜合迭代方法求解函數(shù),調(diào)用后求解;

(2) 使用steffenson法求解。


1、牛頓迭代法


  • 又稱(chēng)為牛頓-拉弗森方法(Newton-Raphson method),單變量下又稱(chēng)為切線法。它是一種在實(shí)數(shù)域和復(fù)數(shù)域上近似求解方程的方法。方法使用函數(shù)f (x)的泰勒級(jí)數(shù)的前面幾項(xiàng)來(lái)尋找方程f (x) = 0的根。用牛頓迭代法解非線性方程,是把非線性方程f(x) = 0線性化的一種近似方法。

  • 牛頓迭代法的本質(zhì)是一種線性化方法,其基本思想是將非線性方程f(x)=0逐步歸結(jié)為某種線性方程來(lái)求解。

    設(shè)方程f(x)=0有近似根,將函數(shù)f(x)在點(diǎn)處展開(kāi),則有

    image15885805357980.png

    于是,方程f(x)=0就可以近似表示為

image15885805567920-1588580565824.png

對(duì)該方程求解,得到
x_{k+1}=x_k-\frac{f(x_k)}{f^{'}(x_k)} \quad\quad (k=0,1,2,3...)
該迭代公式即為 牛頓迭代公式。

x_{k+1}y=f(x)的切線在點(diǎn)(x_k,f(x_k))處與x軸交點(diǎn)的橫坐標(biāo)。

image-20200504153139230.png

算法流程圖:

image-20200504163537448.png

對(duì)牛頓迭代法更加直觀的理解

NewtonIteration_Ani.gif

2、牛頓下山法


在牛頓迭代法中,有時(shí)候會(huì)出現(xiàn)如下迭代回退的情況,導(dǎo)致出現(xiàn)死循環(huán)

20170531114836051.png

當(dāng)選取初值有困難時(shí),可改用如下迭代格式,以擴(kuò)大初值的選取范圍,
x_{n+1}=x_n-\lambda\frac{f(x_k)}{f^{'}(x_k)}
其中λ稱(chēng)為下山因子,λ選取應(yīng)滿足單調(diào)性條件
|f(x_{n+1}|<|f(x_n)|
這樣的算法稱(chēng)下山法.將下山法和牛頓法結(jié)合起來(lái)使用的方法,稱(chēng)為牛頓下山法.

下山因子λ的選擇是逐步探索的過(guò)程.從λ=1開(kāi)始反復(fù)將λ減半進(jìn)行試算,如果能定出λ使單調(diào)性條件成立則“下山成功”.與此相反,如果找不到使單調(diào)性條件成立的λ,則“下山失敗”.此時(shí)需另選初值x0重算

3、割線法引入


割線法又稱(chēng)為弦截法法

  • 為什么要引入割線法?

    在牛頓下山法中我們需要計(jì)算被切割點(diǎn)的倒數(shù),在多次迭代情況下,調(diào)用MatLab自帶的求導(dǎo)函數(shù)會(huì)占用很多計(jì)算資源,所以,我們采用向后有限差分來(lái)近似導(dǎo)數(shù)計(jì)算,什么意思?

    上圖:

    image-20200504152943768.png

    在割線法中,省略了求導(dǎo)步驟,用x_kx_{k-1}兩點(diǎn)連成的直線與x軸的交點(diǎn)近似代替了牛頓法中x_{k+1}的求解

    由此,減少了計(jì)算資源的消耗

  • 需要注意

    弦截法和牛頓迭代法都是線性化方法,牛頓迭代法在計(jì)算x_{n+1}時(shí)只用到前一步的值x_n,而弦截法用到前面兩步的結(jié)果x_nx_{n+1},因此使用割線法必須先給出兩個(gè)初值值x_0,x_{1}.

牛頓下山法和割線法結(jié)合算法如下:

function [x,iter,X]=newton_secant(fun,x0,x1,eps,maxiter)
% 牛頓下山+割線法求解非線性方程的根
% 輸入?yún)?shù):
%      ---fun:迭代函數(shù)
%      ---x0,x1:初始迭代點(diǎn)
%      ---eps:精度要求,默認(rèn)值為1e-6
%      ---maxiter:最大迭代次數(shù),默認(rèn)值為1e4
% 輸出參數(shù):
%      ---x:非線性方程的近似根
%      ---iter:迭代次數(shù)
%      ---X:每一步迭代的結(jié)果
if nargin<3,error('輸入?yún)?shù)至少需要3個(gè)!'),end
if nargin<4|isempty(eps),eps=1e-6;end
if nargin<5|isempty(maxiter),maxiter=1e4;end
f0 = feval(fun,x0); % 計(jì)算x0處的函數(shù)值
f1 = feval(fun,x1); % 計(jì)算x1處的函數(shù)值
k=0;err=1;%k統(tǒng)計(jì)迭代次數(shù),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);  % 牛頓下山迭代,割線法代替求導(dǎo)
        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('迭代次數(shù)超限,迭代失?。?)
end
x=x2;iter=k;X=X';
end

解決第一個(gè)問(wèn)題

  • 每月還款數(shù)額計(jì)算公式如圖:
    cc11728b4710b912a1749141cefdfc03934522d3.png

P:貸款本金

R:月利率

N:還款期數(shù)

月利率 = 年利率/12

  • 所以我們有

180\times\frac{x(1+x)^{240}}{(1+x)^{240}-1}-1.2=0

fun=@(x)180*(x*(1+x)^240)/((1+x)^240-1)-1.2;

命令行調(diào)用函數(shù)解方程:

[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\%

4、steffenson法

Steffeson法的導(dǎo)數(shù)近似方法:
f^{'}(x_k)\approx\frac{f(x_k+f(x_k))-f(x_k)}{f(x_k)}
則迭代公式為:
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...)
廢話少說(shuō),上代碼:

function [x,iter,X] = steffenson(fun,x0,eps,maxiter)
% Steffenson法求解非線性方程的根
% 輸入?yún)?shù):
%      ---fun:迭代函數(shù)
%      ---x0:初始迭代點(diǎn)
%      ---eps:精度要求,默認(rèn)值為1e-6
%      ---maxiter:最大迭代次數(shù),默認(rèn)值為1e4
% 輸出參數(shù):
%      ---x:非線性方程的近似根
%      ---iter:迭代次數(shù)
%      ---X:每一步迭代的結(jié)果
if nargin<2,error('輸入?yún)?shù)至少需要2個(gè)!'),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); % 計(jì)算x0處的函數(shù)值
    x1=x0-f0^2/(feval(fun,x0+f0)-f0);  % Steffenson法迭代公式
    err=x1-x0;
    x0 = x1; % 更新x0數(shù)值
    X(k)=x1;
end
x=x1;iter=k;X=X';

命令行調(diào)用函數(shù)解方程:

[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\%,但是迭代次數(shù)為6次

二、非線性方程組求解

通過(guò)以下問(wèn)題學(xué)習(xí)此知識(shí)點(diǎn):

用牛頓法求解二元方程組的根
\begin{array}{cc} x^2 \mathrm{cos2x}+y^2 \mathrm{sin(2y)}=1 & \\ x^3 +y^3 -6\mathrm{cos(2xy)}=-1 & \end{array}

1、牛頓法解方程組


非線性方程組的求法有很多,此處僅對(duì)使用牛頓法求解非線性方程組的根進(jìn)行學(xué)習(xí)。
將多元向量函數(shù)F(x) 在點(diǎn)處展開(kāi)
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)的Jacobi矩陣

因此,可以得到求解非線性方程組的迭代方程
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)
---------流程圖--------

image-20200504181653878.png

上代碼:

function [x,iter,X]=newtong(fun,x0,eps,maxiter)
% Newton法求解非線性方程組的根
% 輸入?yún)?shù):
%      ---fun:迭代函數(shù)
%      ---x0:初始迭代點(diǎn)向量
%      ---eps:精度要求,默認(rèn)值為1e-6
%      ---maxiter:最大迭代次數(shù),默認(rèn)值為1e4
% 輸出參數(shù):
%      ---x:非線性方程的近似解向量
%      ---iter:迭代次數(shù)
%      ---X:每一步迭代的結(jié)果
if nargin<2,error('輸入?yún)?shù)至少需要2個(gè)!'),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);  % 求函數(shù)fun的值和jacobi矩陣
    x1=x0-J\fx0;  % 牛頓法迭代公式
    err=norm(x1-x0);
    x0=x1;
    X(k,:)=x1;
end
if k==maxiter
    error('迭代次數(shù)超限,迭代失??!')
end
x=x1;iter=k;
end


function [y,J]=fun(x)
    % 非線性方程組
    % 函數(shù)文件描述,返回函數(shù)值和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;  % 聲明符號(hào)變量
    %J=jacobian([2*xx-yy-exp(-xx);-xx+2*yy-exp(-yy)],[xx yy]);  % 求符號(hào)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
          ……
最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點(diǎn),簡(jiǎn)書(shū)系信息發(fā)布平臺(tái),僅提供信息存儲(chǔ)服務(wù)。

友情鏈接更多精彩內(nèi)容