利用牛頓迭代法求解非線性方程組

       近期一個哥們。是用牛頓迭代法求解一個四變量方程組的最優解問題，從網上找了代碼去改進。可是總會有點不如意的地方。迭代的次數過多。可是卻沒有提高精度，真是令人揪心。

       經分析，發現是這個方程組中存在非常多局部的極值點，是用牛頓迭代法不能不免進入局部極值的問題，更程序的初始值有關！

       發現自己好久沒有是用Matlab了。順便從網上查了查代碼，自己來改動一下！

先普及一下牛頓迭代法：(來自百度百科)

       牛頓迭代法（Newton's method）又稱為牛頓-拉夫遜（拉弗森）方法（Newton-Raphson method），它是牛頓在17世紀提出的一種在實數域和復數域上近似求解方程的方法。

多數方程不存在求根公式，因此求精確根很困難，甚至不可能，從而尋找方程的近似根就顯得特別重要。方法使用函數f(x）的泰勒級數的前面幾項來尋找方程f(x) = 0的根。牛頓迭代法是求方程根的重要方法之中的一個。其最大長處是在方程f(x) = 0的單根附近具有平方收斂，並且該法還能夠用來求方程的重根、復根，此時線性收斂，可是可通過一些方法變成超線性收斂。另外該方法廣泛用於計算機編程中。

        設r是f(x)=0的根。選取x0作為r的初始近似值。過點(x0,f(x0))做曲線的切線，求出該切線與x軸的交點，並求出該點的橫坐標，稱作x1是r的一次近似。如此就能夠推導出牛頓迭代公式。

         已經證明。假設是連續的。而且待求的零點是孤立的，那么在零點周圍存在一個區域。僅僅要初始值位於這個鄰近區域內，那么牛頓法必然收斂。 而且，假設不為0, 那么牛頓法將具有平方收斂的性能. 粗略的說。這意味着每迭代一次，牛頓法結果的有效數字將添加一倍。

        在網上查了一些代碼。都是能指定某幾個函數進行求導的。並且要是改變函數的個數，卻又要對原始程序大動干戈。真的是揪心。

        找到了http://hi.baidu.com/aillieo/item/f4d2c4de85af6be954347f25 這個程序。貌似在Matlab上不能非常好的執行，對於數據的返回值為空沒有做處理。后來又找了一個網易朋友的博客，將他的代碼拿過來跑跑，還能夠，可是對於不同的函數方程組，以及變量個數就不同了，真的是揪心。這個就是程序設計和編碼的問題了！

       自己就拿來改了改，能夠支持多方程組和多變量了！

以下附上我的代碼！求大家指導！

 

function [r,n]=mulNewton(x0,funcMat,var,eps)
% x0為兩個變量的起始值,funcMat是兩個方程,var為兩個方程的兩個變量,eps控制精度
% 牛頓迭代法解二元非線性方程組
if nargin==0
    x0 = [0.2,0.6];
    funcMat=[sym('(15*x1+10*x2)-((40-30*x1-10*x2)^2*(15-15*x1))*5e-4')...
             sym('(15*x1+10*x2)-((40-30*x1-10*x2)*(10-10*x2))*4e-2')];
    var=[sym('x1') sym('x2')];
    eps=1.0e-4;
end

n_Var = size(var,2);%變量的個數
n_Func = size(funcMat,2);%函數的個數
n_X = size(x0,2);%變量的個數

if n_X ~= n_Var && n_X ~= n_Func
    fprintf('Expression Error!\n');
    exit(0);
end

r=x0-myf(x0, funcMat, var)*inv(dmyf(x0, funcMat, var));
n=0;
tol=1;
while tol>=eps
    x0=r;
    r=x0-myf(x0, funcMat, var)*inv(dmyf(x0, funcMat, var));
    tol=norm(r-x0);
    n=n+1;
    if(n>100000)
        disp('迭代步數太多，方程可能不收斂');
        return;
    end
end
end % end mulNewton

 

function f=myf(x,funcMat, varMat)
% 輸入參數x為兩個數值,func為1*2符號變量矩陣,var為1*2符號變量矩陣中的變量
% 返回值為1*2矩陣,內容為數值

n_X = size(x,2);%變量的個數
f_Val = zeros(1,n_X);
for i=1:n_X
    tmp_Var = cell(1,n_X);
    tmp_X = cell(1,n_X);
    for j=1:n_X
        tmp_Var{j} = varMat(1,j);
        tmp_X{j} = x(1,j);
    end
    f_Val(i) = subs(funcMat(1, i), tmp_Var, tmp_X);
end
f=f_Val;
end % end myf

function df_val=dmyf(x, funcMat, varMat)
% 返回值為2*2矩陣,內容為數值
%df=[df1/x1, df1/x2;
%    df2/x1. df2/x2];
n_X = size(x,2);%變量的個數
df =cell(n_X, n_X);
for i=1:n_X
    for j=1:n_X
        df{i,j} = diff(funcMat(1, i), varMat(1, j));
    end
end

df_val=zeros(n_X, n_X);

for i=1:n_X
    for j=1:n_X
        tmp_Var = cell(1,n_X);
        tmp_X = cell(1,n_X);
        for k=1:n_X
            tmp_Var{k} = varMat(1,k);
            tmp_X{k} = x(1,k);
        end
        df_val(i,j) = subs(df{i,j}, tmp_Var, tmp_X);
    end
end
end % end dmyf