牛顿迭代法，又名切线法，这里不详细介绍，简单说明每一次牛顿迭代的运算：首先将各个方程式在一个根的估计值处线性化（泰勒展开式忽略高阶余项），然后求解线性化后的方程组，最后再更新根的估计值。下面以求解最简单的非线性二元方程组为例（平面二维定位最基本原理），贴出源代码：
1、新建函数fun.m，定义方程组
1 function f=fun(x); 2 %定义非线性方程组如下 3 %变量x1 x2 4 %函数f1 f2 5 syms x1 x2 6 f1 = sqrt((x1-4)^2 + x2^2)-sqrt(17); 7 f2 = sqrt(x1^2 + (x2-4)^2)-5; 8 f=[f1 f2];
2、新建dfun.m，求出一阶微分方程
1 function df=dfun(x); 2 f=fun(x); 3 df=[diff(f,'x1');diff(f,'x2')]; %雅克比矩阵
3、建立newton.m，执行牛顿迭代过程
1 clear;clc 2 format; 3 x0=[0 0]; % 迭代初始值 4 eps = 0.00001; % 定位精度要求 5 for i = 1:10 6 f = double(subs(fun(x0),{'x1' 'x2'},{x0(1) x0(2)})); 7 df = double(subs(dfun(x0),{'x1' 'x2'},{x0(1) x0(2)})); % 得到雅克比矩阵 8 x = x0 - f/df; 9 if(abs(x-x0) < eps) 10 break; 11 end 12 x0 = x; % 更新迭代结果 13 end 14 disp('定位坐标：'); 15 x 16 disp('迭代次数：'); 17 i
结果如下：
定位坐标：
x =
0.0000 -1.0000
迭代次数：
i =
4