MATLAB求解非线性方程组

function y=f_sb(x) %构造线性方程组
p=1025;
g=9.8;
y(1)=-p*g*x(1)*pi+12000+sin(x(2))*x(3);
y(2)=-x(3)*cos(x(2))+0.625*pi*(2-x(1))*x(4)*x(4);

y(3)=-x(3)*sin(x(2))-pi*0.025*0.025*p*g+100+x(5)*sin(x(6));
y(4)=-x(3)*cos(x(2))+x(5)*cos(x(6));

y(5)=-x(5)*sin(x(6))-pi*0.025*0.025*p*g+100+x(7)*sin(x(8));
y(6)=-x(5)*cos(x(6))+x(7)*cos(x(8));

y(7)=-x(7)*sin(x(8))-pi*0.025*0.025*p*g+100+x(9)*sin(x(10));
y(8)=-x(7)*cos(x(8))+x(9)*cos(x(10));

y(9)=-x(9)*sin(x(10))-pi*0.025*0.025*p*g+100+x(11)*sin(x(12));
y(10)=-x(9)*cos(x(10))+x(11)*cos(x(12));
a=0.625*x(4)*x(4)*pi*(2-x(1));
y(11)=-0.15*0.15*p*g-x(11)*sin(x(12))+1000+0.625*x(4)*x(4)*pi*(2-x(1))*cosh(x(13)/a)*sin(x(14));
y(12)=-x(11)*cos(x(12))+0.625*x(4)*x(4)*pi*(2-x(1))*cosh(x(13)/a)*cos(x(14));
y(13)=-a*sinh(x(13)/a)+22.05;
y(14)=-x(1)-(sin(x(2))+sin(x(6))+sin(x(8))+sin(x(10))+sin(x(12)))-x(15)+18;
y(15)=-x(15)+a-cosh(x(13)/a-1);
y(16)=-p*(pi*x(1)+0.025*0.025*pi*4+0.15*0.15*pi)+1000+40+100+1200;    

options=optimset('Display','off');
x=fsolve(@f_sb,ones(15,1)',options) %求解
y=f_sb(x) %验证

这是最基本的方法，需要给定初值，误差较大。