求助matlab与1stopt求解非线性方程组

wangbo1113

我是一名博士研究生，一直用matlab编程求解问题，最近遇到一个非线性方程组迭代求解问题，由于初值无法确定，所以困扰了好长时间；在论坛中看到可以用1stopt求解非线性方程组无需给定初值，所以也尝试了一下，但是由于是新手，所以对于程序中的许多简单问题都无法确定。所以向各位大侠求教。
下面将我的问题描述一下：在转速n为循环变量的前提下，给定a1，a2，a3一组初值，求解一个4元非线性方程组得出X1、X2、X3、X4，（其中，并且每个X都是一个16维数组）；利用求出的X（4*16）带入另外3个方程组中（每个方程中有求和），重新求解a1，a2，a3，直到满足精度要求。
针对上述问题，我编写了matlab程序，由于初值问题，总是算不出貌似合理的结果。我把程序贴出来，请大侠指教。
%% 主函数
begin=3000;  step=3000;   finish=3000;% 变参表示转速
csstart=1;  csstep=1;csend=16;
XX=zeros(4,length(csstart:csstep:csend));
a=[0.02,0.0,0/180*pi];%初始delta_a,delta_r,theta
x0=zeros(4,16);x0(1:2,:)=0.4;n=0;
for biancan=begin:step:finish
    m=0;
    for csi=csstart:csstep:csend
        m=m+1;
        options = optimset('MaxFunEvals',10000000,'MaxIter',10000);
        x = fsolve(@(x) highspeed_displacement_4(x,biancan,csi,a),[0.7,0.2,0.1,0.1],options);
         XX(:,m)=x';
    end
    y = fsolve(@(y) highspeed_later_3j(biancan,XX,y),a,options);
   
    while(abs(y(1)-a(1))>1e-2)
        n=n+1
        a(1)=y(1);m=0;
        for csi=csstart:csstep:csend
            m=m+1;
            x = fsolve(@(x) highspeed_displacement_4(x,biancan,csi,a),[0.4,0.4,0,0],options);
            XX(:,m)=x';
        end
         y = fsolve(@(y) highspeed_later_3j(biancan,XX,y),a,options);
    end
end
function F=highspeed_displacement_4(x,biancan,csi,a)
D=23;            Z=16;       alpha0=40/180*pi;   dm=125;                             rho=7850;          m=rho*pi*D^3/6;     J=rho*pi*D^5/60;
fi=0.515;        fo=0.525;   B=fi+fo-1;          Re=0.5*dm+(fi-0.5)*D*cos(alpha0);   gamma_pie=D/dm;
lambda_ij=0;     lambda_oj =2;%外滚道控制
% lambda_ij=1;     lambda_oj =1;%内外滚道平均分担
speed_n=biancan;       omega=2*pi*speed_n/60;
j=csi;
psi_j=2*pi*(j-1)/Z;
A1j=B*D*sin(alpha0)+a(1)+Re*a(3)*cos(psi_j);
A2j=B*D*cos(alpha0)+a(2)*cos(psi_j);
%%  表示cosαoj，sinαoj，cosαij，sinαij
COS_alpha_oj=x(2)/((fo-0.5)*D+x(4));
SIN_alpha_oj=x(1)/((fo-0.5)*D+x(4));
COS_alpha_ij=(A2j-x(2))/((fi-0.5)*D+x(3));
SIN_alpha_ij=(A1j-x(1))/((fi-0.5)*D+x(3));
%% 求载荷位移系数K  Q=K*delta^1.5
%   钢球与内圈接触————————*————————————*————————（%钢球与外圈接触 %r_x2=(dm/cos(alpha/180*pi)+D)/2;    r_y2=D*fo;）
% Koj=Q_delta(D,fi,fo,dm,x(2),1);   Kij=Q_delta(D,fi,fo,dm,x(1),1);
% 求Kij
r_x1=D/2;    r_y1=D/2;%钢球半径
r_x2=(dm/COS_alpha_ij-D)/2;    r_y2=D*fi;
rho_x1=1/r_x1;rho_x2=1/r_x2;rho_y1=1/r_y1;rho_y2=1/r_y2;
rho_zong=rho_x1+rho_y1+rho_x2+rho_y2;
Rx=1/(rho_x1+rho_x2);  Ry=1/(rho_y1+rho_y2);
kappa=1.0339*(Ry/Rx)^0.636;  EE=1.0003+0.5968*Rx/Ry;  FF=1.5277+0.6023*log(Ry/Rx);
delta_8=2*FF/pi*(pi/(2*kappa^2*EE))^(1/3);
Kij=2.15*10^5*rho_zong^(-0.5)*delta_8^(-1.5);
%求Koj
r_x2=(dm/COS_alpha_oj-D)/2;
rho_x1=1/r_x1;rho_x2=1/r_x2;rho_y1=1/r_y1;rho_y2=1/r_y2;
rho_zong=rho_x1+rho_y1+rho_x2+rho_y2;
Rx=1/(rho_x1+rho_x2);  Ry=1/(rho_y1+rho_y2);
kappa=1.0339*(Ry/Rx)^0.636;  EE=1.0003+0.5968*Rx/Ry;  FF=1.5277+0.6023*log(Ry/Rx);
delta_8=2*FF/pi*(pi/(2*kappa^2*EE))^(1/3);
Koj=2.15*10^5*rho_zong^(-0.5)*delta_8^(-1.5);
%%
omega_m_ratio=(1-gamma_pie*COS_alpha_ij)/(1+COS_alpha_ij*COS_alpha_oj+SIN_alpha_oj*SIN_alpha_ij);
beta=atan(SIN_alpha_oj/COS_alpha_oj+gamma_pie);
omega_r_ratio=-1/((COS_alpha_oj+tan(beta)*SIN_alpha_oj)/(1+gamma_pie*COS_alpha_oj)+((COS_alpha_ij+tan(beta)*SIN_alpha_ij)/(1-gamma_pie*COS_alpha_ij))*gamma_pie*cos(beta));%内滚道旋转
% omega_r_ratio=1/(cos(x(2))+tan(beta)*sin(x(2)))/(1+gamma_pie*cos(x(2)))+((cos(x(1))+tan(beta)*sin(x(1)))/(1-gamma_pie*cos(x(1))))/gamma_pie/cos(beta);%外滚道旋转
Mgj=J*omega_r_ratio*omega_m_ratio*omega^2*sin(beta);% 陀螺力矩
Fcj=0.5*m*dm*omega^2*omega_m_ratio^2; %离心力
%%  函数
F1=(A1j-x(1))^2+(A2j-x(2))^2-((fi-0.5)*D+x(3))^2;
F2=x(1)^2+x(2)^2-((fo-0.5)*D+x(4))^2;
F3=(lambda_oj*Mgj*x(2)/D-Koj*x(4)^1.5*x(1))/((fo-0.5)*D+x(4))+(Kij*x(3)^1.5*(A1j-x(1))-lambda_ij*Mgj*(A2j-x(2))/D)/((fi-0.5)*D+x(3));
F4=(lambda_oj*Mgj*x(1)/D+Koj*x(4)^1.5*x(2))/((fo-0.5)*D+x(4))-(Kij*x(3)^1.5*(A2j-x(2))+lambda_ij*Mgj*(A1j-x(1))/D)/((fi-0.5)*D+x(3))-Fcj;
F=[F1 F2 F3 F4 ]';
function F=highspeed_later_3j(biancan,a,y)
%% 变量
% x为highspeed_f_4计算出来的X1j，X2j，delta_ij ,delta_oj
% biancan 仍然为转速 r/min
% x=XX(3:end,:);
% syms delta_a delta_r  theta
%% 参数输入
D=23;            Z=16;       alpha0=40/180*pi;   dm=125;                             rho=7850;          m=rho*pi*D^3/6;     J=rho*pi*D^5/60;
fi=0.515;        fo=0.525;   B=fi+fo-1;          Re=0.5*dm+(fi-0.5)*D*cos(alpha0);   gamma_pie=D/dm;
lambda_ij=0;     lambda_oj =2;%外滚道控制
% lambda_ij=1;     lambda_oj =1;%内外滚道平均分担
speed_n=biancan;       omega=2*pi*speed_n/60;
Fa=10000;Fr=0;M=0;F5=0;F6=0;F7=0;
%%
for j=1:Z
    psi_j=2*pi*(j-1)/Z;
    A1j=B*D*sin(alpha0)+y(1)+Re*y(3)*cos(psi_j);
    A2j=B*D*cos(alpha0)+y(2)*cos(psi_j);
    %%  表示cosαoj，sinαoj，cosαij，sinαij
    COS_alpha_oj=a(2,j)/((fo-0.5)*D+a(4,j));
    SIN_alpha_oj=a(1,j)/((fo-0.5)*D+a(4,j));
    COS_alpha_ij=(A2j-a(2,j))/((fi-0.5)*D+a(3,j));
    SIN_alpha_ij=(A1j-a(1,j))/((fi-0.5)*D+a(3,j));
    %% 求载荷位移系数K  Q=K*delta^1.5
    r_x1=D/2;    r_y1=D/2;%钢球半径
    r_x2=(dm/COS_alpha_ij-D)/2;    r_y2=D*fi;
    rho_x1=1/r_x1;rho_x2=1/r_x2;rho_y1=1/r_y1;rho_y2=1/r_y2;
    rho_zong=rho_x1+rho_y1+rho_x2+rho_y2;
    Rx=1/(rho_x1+rho_x2);  Ry=1/(rho_y1+rho_y2);
    kappa=1.0339*(Ry/Rx)^0.636;  EE=1.0003+0.5968*Rx/Ry;  FF=1.5277+0.6023*log(Ry/Rx);
    delta_8=2*FF/pi*(pi/(2*kappa^2*EE))^(1/3);
    Kij=2.15*10^5*rho_zong^(-0.5)*delta_8^(-1.5);
    %%
    omega_m_ratio=(1-gamma_pie*COS_alpha_ij)/(1+COS_alpha_ij*COS_alpha_oj+SIN_alpha_oj*SIN_alpha_ij);
    beta=atan(SIN_alpha_oj/COS_alpha_oj+gamma_pie);
    omega_r_ratio=-1/((COS_alpha_oj+tan(beta)*SIN_alpha_oj)/(1+gamma_pie*COS_alpha_oj)+((COS_alpha_ij+tan(beta)*SIN_alpha_ij)/(1-gamma_pie*COS_alpha_ij))*gamma_pie*cos(beta));%内滚道旋转
    % omega_r_ratio=1/(cos(x(2))+tan(beta)*sin(x(2)))/(1+gamma_pie*cos(x(2)))+((cos(x(1))+tan(beta)*sin(x(1)))/(1-gamma_pie*cos(x(1))))/gamma_pie/cos(beta);%外滚道旋转
    Mgj=J*omega_r_ratio*omega_m_ratio*omega^2*sin(beta);% 陀螺力矩
    %%  函数
    f5=(Kij*(A1j-a(1,j))*a(3,j)^1.5-lambda_ij*Mgj*(A2j-a(2,j))/D)/((fi-0.5)*D+a(3,j));                                  F5=F5+f5;
    f6=(Kij*(A2j-a(2,j))*a(3,j)^1.5-lambda_ij*Mgj*(A1j-a(1,j))/D)/((fi-0.5)*D+a(3,j))*cos(psi_j);                       F6=F6+f6;
    f7=((Kij*(A1j-a(1,j))*a(3,j)^1.5-lambda_ij*Mgj*(A2j-a(2,j))/D)*Re/((fi-0.5)*D+a(3,j))+lambda_ij*fi*Mgj)*cos(psi_j); F7=F7+f7;
end
F5=Fa-F5;    F6=Fr-F6;        F7=M-F7;
F=[F5,F6,F7]';
然后我也编了求解X1-X4的1stopt程序，但是总有一些过程变量，不知道带入值了没有，虽然算出了结果但是不跟确实是否正确，也求助。
Parameter  x1,x2,x3,x4;
Constant  D=23;               Constant Z=16;              Constant dm=125;           Constant  alpha0=40/180*pi;
Constant  rho=7850;           Constant m=rho*pi*D^3/6;    Constant JJ=rho*pi*D^5/60;
Constant  fi=0.515;           Constant fo=0.525;          Constant  lambda_ij=0;     Constant lambda_oj =2;
Constant  B=fi+fo-1;          Constant Re=0.5*dm+(fi-0.5)*D*cos(40/180*pi);             Constant gamma_pie=D/dm;
Constant  speed_n=0;                                      Constant j=1;
Constant omega=2*pi*speed_n/60,                           Constant psi=2*pi*(j-1)/Z;
Constant a1=0.025;            Constant a2=0;              Constant a3=0*pi/180;
Constant AA1=B*D*sin(40/180*pi)+a1+Re*a3*Cos(2*pi*(j-1)/Z) ;   //
Constant AA2=B*D*cos(40/180*pi)+a2*cos(2*pi*(j-1)/Z);  //
ConstStr COS_alpha_oj=x2/((fo-0.5)*D+x4);
ConstStr SIN_alpha_oj=x1/((fo-0.5)*D+x4);
ConstStr COS_alpha_ij=(AA2-x2)/((fi-0.5)*D+x3);
ConstStr SIN_alpha_ij=(AA1-x1)/((fi-0.5)*D+x3);
Constant r_x1=D/2;         Constant  r_y1=D/2;          Constant    r_y2=D*fi;
Constant rhox1=2/D;        Constant  rhoy1=2/D;         Constant    rhoy2=1/(D*fi);    Constant   Ry=1/(2/D+1/(D*fi));
ConstStr r_x2=(dm/COS_alpha_ij-D)/2,
         rhox2=2/(dm/COS_alpha_ij-D),
         rhozong=rhox1+rhoy1+rhox2+rhoy2,
         Rx=1/(rhox1+rhox2),
         kappa=1.0339*(Ry/Rx)^0.636,
         EE=1.0003+0.5968*Rx/Ry,
         FF=1.5277+0.6023*log(Ry/Rx),
         delta_8=2*FF/pi*(pi/(2*kappa^2*EE))^(1/3),
         Kij=2.15*10^5*rhozong^(-0.5)*delta_8^(-1.5);
ConstStr r_x22=(dm/COS_alpha_oj-D)/2,
         rhox22=2/(dm/COS_alpha_oj-D),
         rhozong2=rhox1+rhoy1+rhox22+rhoy2,
         Rx2=1/(rhox1+rhox22),
         kappa2=1.0339*(Ry/Rx2)^0.636,
         EE2=1.0003+0.5968*Rx2/Ry,
         FF2=1.5277+0.6023*log(Ry/Rx2),
         delta_82=2*FF2/pi*(pi/(2*kappa2^2*EE2))^(1/3),
         Koj=2.15*10^5*rhozong2^(-0.5)*delta_82^(-1.5);
ConstStr omegamratio=(1-gamma_pie*COS_alpha_ij)/(1+COS_alpha_ij*COS_alpha_oj+SIN_alpha_oj*SIN_alpha_ij);
ConstStr beta=atan((x1/((fo-0.5)*D+x4))/(x2/((fo-0.5)*D+x4))+gamma_pie);
ConstStr omegarratio=-1/((COS_alpha_oj+tan(beta)*SIN_alpha_oj)/(1+gamma_pie*COS_alpha_oj)+((COS_alpha_ij+tan(beta)*SIN_alpha_ij)/(1-gamma_pie*COS_alpha_ij))*gamma_pie*cos(beta));
ConstStr Mgj=JJ*omegarratio*omegamratio*omega^2*sin(beta);
ConstStr Fcj=0.5*m*dm*omega^2*omegamratio^2;
Function
(AA1-x1)^2+(AA2-x2)^2-((fi-0.5)*D+x3)^2=0;
x1^2+x2^2-((fo-0.5)*D+x4)^2=0;
(lambda_oj*Mgj*x2/D-Koj*x4^1.5*x1)/((fo-0.5)*D+x4)+(Kij*x3^1.5*(AA1-x1)-lambda_ij*Mgj*(AA2-x2)/D)/((fi-0.5)*D+x3)=0;
(lambda_oj*Mgj*x1/D+Koj*x4^1.5*x2)/((fo-0.5)*D+x4)-(Kij*x3^1.5*(AA2-x2)+lambda_ij*Mgj*(AA1-x1)/D)/((fi-0.5)*D+x3)-Fcj=0;

tkcc

我觉得我们要解的是同一组方程，我也在搞这方面的东西，也困扰在这些方程组上面。

tkcc

不知道你现在做的怎么样了，有时间可以交流一下吗？我的qq 116573877