利用牛顿迭代法求解非线性方程组

       最近一个哥们，是用牛顿迭代法求解一个四变量方程组的最优解问题，从网上找了代码去改进，但是总会有点不如意的地方，迭代的次数过多，但是却没有提高精度，真是令人揪心！

       经分析，发现是这个方程组中存在很多局部的极值点，是用牛顿迭代法不能不免进入局部极值的问题，更程序的初始值有关！

       发现自己好久没有是用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