求助matlab与1stopt求解非线性方程组
我是一名博士研究生，一直用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=;%初始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),,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),,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));
%% 求载荷位移系数KQ=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=';
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_rtheta
%% 参数输入
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));
%% 求载荷位移系数KQ=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=';
然后我也编了求解X1-X4的1stopt程序，但是总有一些过程变量，不知道带入值了没有，虽然算出了结果但是不跟确实是否正确，也求助。
Parameterx1,x2,x3,x4;
ConstantD=23; Constant Z=16; Constant dm=125; Constantalpha0=40/180*pi;
Constantrho=7850; Constant m=rho*pi*D^3/6; Constant JJ=rho*pi*D^5/60;
Constantfi=0.515; Constant fo=0.525; Constantlambda_ij=0; Constant lambda_oj =2;
ConstantB=fi+fo-1; Constant Re=0.5*dm+(fi-0.5)*D*cos(40/180*pi); Constant gamma_pie=D/dm;
Constantspeed_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; Constantr_y1=D/2; Constant r_y2=D*fi;
Constant rhox1=2/D; Constantrhoy1=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;
我觉得我们要解的是同一组方程，我也在搞这方面的东西，也困扰在这些方程组上面。 不知道你现在做的怎么样了，有时间可以交流一下吗？我的qq 116573877
页:
[1]