x0=[0.8, 0.8, 0.2, 0.2]';
 m0=[0.5,0.5]'; 
tol=1e-6;
syms x y z p m1 m2
h1=y-(x^3)-(z^2);
h2=(x^2)-y-p^2;
h=[h1;h2];
f=-x; 
L=f+ m1*h1+m2*h2;
h0 = subs(h, {x,y,z,p}, [x0(1,1), x0(2,1), x0(3,1), x0(4,1)]); 
h0=double(h0);
%the system h evaluated at the intial vector
g1=gradient(h1, [x, y, z, p]);
g2=gradient(h2, [x, y, z, p]);
J=[g1 g2];
J0=subs(J, x, x0(1,1)); 
J0=subs(J0, y, x0(2,1));
J0=subs(J0, z, x0(3,1));
J0=subs(J0, p, x0(4,1));
%J evaluated at the intial vector
J0=double(J0)
n=size(J0,1); %number of rows in J0
m=size(J0,2); %number of columns in J0
DL=gradient(L, [x,y,z,p,m1,m2]); 
DL=subs(DL, y, x0(2,1));
DL=subs(DL, z, x0(3,1));
DL=subs(DL, p, x0(4,1));
DL=subs(DL, x, x0(1,1));
DL=subs(DL,m1, m0(1,1));
DL=subs(DL,m2, m0(2,1));
DL=double(DL)           
DLd=DL(1:n,1);
H=hessian(L, [x,y,z,p]); 
H0 =subs(H,x,x0(1,1));
H0=subs(H0,y,x0(2,1));
H0=subs(H0,z, x0(3,1));
H0=subs(H0,p, x0(4,1));
H0=subs(H0,m1, m0(1,1));
H0=subs(H0,m2, m0(2,1));
H0=double(H0)
[Q, R]= qr(J0);
Y=Q(1:n,1:m);
Z=Q(1:n, m+1:n);
qz=Z'*DLd; qh=h0;
while norm(qz)+ norm(qh) > tol
E=[Z'*H0;J0'];
V=[Z'*DLd;h0];
s0=E\-V
x0=x0+s0 %the new point
J0=subs(J, x, x0(1,1)); 
J0=subs(J0, y, x0(2,1));
J0=subs(J0, z, x0(3,1));
J0=subs(J0, p, x0(4,1));
J0=double(J0);
[Q, R]= qr(J0);
Y=Q(1:n,1:m);
Z=Q(1:n, m+1:n);

r=R(1:m,1:m);
T0=[-1 0 0 0]';
%T evaluated at the new vector x0
m0=r\-(Y'*T0); 
qz=Z'*DLd;
h0 = subs(h, {x,y,z,p}, [x0(1,1), x0(2,1), x0(3,1), x0(4,1)]); 
h0=double(h0);
qh=h0;
end
newvector=[double(x0); double(m0)];

编辑：在代码中进行了一些更正，并在循环内更新了 qz 和 qh。但是，代码似乎永远循环。永远不会违反停止标准，并且 x0 会继续变化，但与 x=(1,1,0,0) 的实际解不同。我将 tol 更改为一个相对较大的值，即 0.7，并且代码生成了一个接近解决方案的输出。当 tol=1e-6 时，双精度会导致这个问题吗？