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

       經分析，發現是這個方程組中存在很多區域性的極值點，是用牛頓迭代法不能不免進入區域性極值的問題，更程式的初始值有關！

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

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

       牛頓迭代法（Newton's method）又稱為牛頓-拉夫遜（拉弗森）方法（Newton-Raphson method），它是牛頓在17世紀提出的一種在

        設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