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