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

我们得首先了解牛顿迭代法的原理是什么：
在这里，我以二阶非线性方程组为例：
f1(x,y)=0
f2(x,y)=0,求解x,y

原理
假设方程组的一组近似解为（x0，y0），将方程f1(x,y)=0与f2(x,y)=0在（x0，y0）处利用二元泰勒级数展开，并取到x，y的一阶近似，则可得到线性方程组：

这两个线性方程组可以写成线性代数里面的矩阵相乘形式：

若系数矩阵的行列式J0不等于0：

最后我们可得到迭代过程中的第一组解：

但是这一组解并不是最优解，我们需要设计一个循环，将这一组得到的x1与y1成为就像第一次循环中的x0和y0，直到得到相邻两个解的差值精度（收敛精度）达到我们预想的值，迭代思路如下：

下面直接上代码，我以解方程

为例

clear all
clc
tol=10^(-8);
x=1;y=1;%初始值
J0=2*x*(x+1)-2*y*(y-3);%求J0
x1=x-((x^2+y^2-5)*(x+1)-((x+1)*y-(3*x+1))*2*y)/J0;%赋值
y1=y-(2*x*((x+1)*y-(3*x+1))-(x^2+y^2-5)*(y-3))/J0;%赋值
while abs(sqrt((x1-x)^2+(y1-y)^2))>tol %开始循环
    x=x1;y=y1;
    J0=2*x*(x+1)-2*y*(y-3);
    x1=x-((x^2+y^2-5)*(x+1)-((x+1)*y-(3*x+1))*2*y)/J0;
    y1=y-(2*x*((x+1)*y-(3*x+1))-(x^2+y^2-5)*(y-3))/J0;
end
x
y