matlab多变量牛顿方法求解非线性方程组

 这本是个作业，但因为各种原因及出错，这个简单的代码写了2个多小时，内心不爽。因此写下博客记录。

clear;
clc;
syms u v
f1 = 6 * u^3 + u * v - 3 * v^3 - 4;
f2 = u^2- 18 * u * v^2 + 16 * v^3 + 1 ;
f=[f1 f2 ];
df0=[diff(f,u);diff(f,v)];
df = df0.';
x0=[-1 -1]; %%%As shown in the book, the solution we find by Newton's way is based on the first estimate
N=200;
% 求解
for i=1:N; % 运算次数
    p=subs(f,{u ,v },{x0(1) x0(2)});% 为subs函数的赋值运算的矩阵
    q=subs(df,{u ,v },{x0(1) x0(2) });% 为subs函数的赋值运算 2*2的矩阵
    x =x0 - grass(q,p);% 得到的一组新数值矩阵 理解为 新的x0
   if norm(x-x0)<eps
      break;
   end % 条件判断 得出的新值是否满足精确度的要求 可以求x矩阵元素平方和的开根
   x0 = vpa(x)% 当条件满足 输出求得的x0
end
x = vpa(x0)

上述的N其实可以根据需要调整，而且初始估计也是对于不同的问题而不同的

其中调用的grass（）函数为朴素的高斯消元法：
function xc = grass(a,b)
n = 2;
for j = 1 : n-1
    if abs(a(j,j))<eps;
    error('zero pivot encountered');
    end
    for i = j+1 : n
        mult = a(i,j)/a(j,j);
        for k = j+1:n
            a(i,k) = a(i,k) - mult*a(j,k);
        end
        b(i) = b(i) - mult*b(j);
    end
end
for i = n : -1 : 1
    for j = i+1 : n
        b(i) = b(i) - a(i,j)*x(j);
    end
    x(i) = b(i)/a(i,i);
end
xc = x;
end