matlab牛顿迭代求解，计算出来复数是怎么回事？

ME!

ME!

x0=x0-s   
            eps=norm(s);
          end
         
    ae=x0(1)
    ai=x0(2)

tammy

把代码贴全，现在给出的代码找不出来问题所在

ME!

这是所有的程序

ME!

可能看不太清楚

朋友

我是来学习的

ME!

yanyongju

可能是数值计算误差，导致刚度矩阵不对称，进而出现复值，建议先检查刚度矩阵

ME!

刚度矩阵时对称的，我在子函数里面看了，我换了几个轴承的尺寸参数，就没有出现复数了，这是怎么回事？

ME!

还有那个迭代的初始值怎么确定啊，好像有时候和这个也有关系

ChaChing

个人水平/时间有限, 楼主的程式对我而言真的觉得有点乱
真的建议看下精华帖,  如何聪明的发帖!

deltai1=sQi^(2/3)*ski

此式当sQi为负数时, deltai1将为复数
至於sQi该不该为负数, 楼主应该比看官更清楚才对

ME!

新手，不懂什么叫精华帖，哪里可以看到这些，
这样是不是清楚一些呢
function Z=dingyayujin7003(n,f)
fa=300;
n=10000;
z=20;%%%%%%%%%%%%%%%%%%滚球数目
dm=(17+35)/1000;%%%轴承中径
db=5.4/1000;%%%%滚球直径
den=3.2e3; %%%%%%%%%%%%%%%%%%%%钢球密度
v1=0.3;%%%%%%%%%%%%%钢的泊松比
v2=0.26;
E1=2.06e11;%%%%%%%%%%钢的弹性模量
E2=3.2e11;
DE=(1-v1^2)/E1+(1-v2^2)/E2;
a0=15*pi/180;%%%初始接触角
rr=db*cos(a0)/dm;%%%%%%%计算主曲率和的所需的中间变量
ffi=0.515;%%%%%%%%%%%%%%%%%%%%内沟道半径系数
ffe=0.525;%%%%%%%%%%%%%%%%%%%%外沟道半径系数(一般外沟道曲率半径要大一些)
ri=ffi*db;   %%%%%%%%%内沟道半径
re=ffe*db;   %%%%%%%%%外沟道半径
PI=(4-(1/ffi)+2*rr/(1-rr))/db;       %%%%钢球和内圈主曲率和p11+p12+p21+p22      p表示曲率半径
FPI=(1/ffi+2*rr/(1-rr))/(4-(1/ffi)+2*rr/(1-rr));  %钢球和内圈主曲率函数|(p11-p12)+(p21-p22)|/(p11+p12+p21+p22)
PE=(4-1/ffe-2*rr/(1+rr))/db;  %%钢球和外圈主曲率和
FPE=(1/ffe-2*rr/(1+rr))/(4-1/ffe-2*rr/(1+rr)); %%钢球和外圈主曲率函数

r=db/dm;
I=(1/60)*den*pi*db^5; %滚动体的转动惯量
Ri=dm/2+(ffi-0.5)*db*cos(a0);%%%%内圈沟道曲率中心轨迹的半径
Rye=(ffe*db)/(2*ffe-1);%%%%计算椭圆积分需要的中间变量
Ryi=(ffi*db)/(2*ffi-1);%%%%计算椭圆积分需要的中间变量

     
     w=n*2*pi/60;      %%%sae,sai的迭代初值越接近最终的值越好
     sae=1*pi/180;     %带s的为迭代中间变量，外接触角ae;
     sai=70*pi/180;   %迭代中间变量ai
     x0=[sae sai]';
          eps=5;%%%%%%%%%标志迭代误差的变量
          while abs(eps)>1e-5 %%%%%%%%%%%%%%%迭代算法
            sae=x0(1);   
            sai=x0(2);         
            
            sRxe=(db*(dm+db*cos(sae)))/(2*dm);%%%%%%椭圆积分的拟合参数
            sRxi=(db*(dm-db*cos(sai)))/(2*dm);%%%%%%椭圆积分的拟合参数
            
            seli1e=1.5277+0.6023*log(Rye/sRxe);%第一类椭圆积分近似计算公式
            seli1i=1.5277+0.6023*log(Ryi/sRxi); %第一类椭圆积分近似计算公式   
            seli2e=1.0003+0.5968*(sRxe/Rye);    %第二类椭圆积分
            seli2i=1.0003+0.5968*(sRxi/Ryi); %第二类椭圆积分     
     
               
            selie=1.0339*(Rye/sRxe)^0.636;    %外圈椭圆率参数，为接触椭圆长半轴与短半轴之比
            selii=1.0339*(Ryi/sRxi)^0.636;  %椭圆率参数，为接触椭圆长半轴与短半轴之比   
            
            
            
            
            sb=atan(sin(sae)/(cos(sae)+r));      %滚动体姿态角（自转轴线与x y坐标轴夹角）对应袁卫公式3.27  
            swb=-w/(((cos(sae)+tan(sb)*sin(sae))/(1+r*cos(sae))+(cos(sai)+tan(sb)*sin(sai))/(1-r*cos(sai)))*r*cos(sb));%滚动体的自转角速度  
            swm=w*(1-r*cos(sai))/(1+cos(sai-sae));%滚动体的公转角速度     
            sfc=(pi/12)*den*dm*db^3*swm^2*w; %滚动体离心力   对应纪宗辉2.7
            smg=I*sin(sb)*swb*swm; %考虑滚动体相对于内套圈的陀螺力矩   
            
   
            ssre1=db*cos(sae)/dm;    %中间变量
            ssri1=db*cos(sai)/dm;    %中间变量
            sPe=(4-(1/ffe)+(2*ssre1)/(1-ssre1))/db;     %滚动体与内套圈接触的曲率和     
            sPi=(4-(1/ffi)+(2*ssri1)/(1-ssri1))/db;      %滚动体与外套圈接触的曲率和
            
            ske=seli1e*((9/2/seli2e/sPe)*(1/pi/selie)^2)^(1/3); %滚动体和外套圈的接触刚度系数 对应公式2.24
            ski=seli1i*((9/2/seli2i/sPi)*(1/pi/selii)^2)^(1/3);   %滚动体和内套圈的接触刚度系数  对应公式2.24
            sQi=fa/(z*sin(sai));    %滚动体和内圈的接触力
            sQe=(fa/z+2*smg*cos(sae)/db)/sin(sae);    %滚动体和外圈的接触力
            
            deltai1=sQi^(2/3)*ski;     %滚动体和内圈接触的变形与载荷的关系，对应公式2.28*P13
            deltae1=sQe^(2/3)*ske;     %滚动体和内圈接触的变形与载荷的关系，对应公式2.29*P13
            F1='(re-0.5*db+deltae1)*cos(sae)+(ri-0.5*db+deltai1)*cos(sai)-(ri+re-db)*cos(a0)'%滚动体和内圈接触的变形与载荷的关系，对应公式2.21*P13
            F2='1/tan(sai)-1/tan(sae)+(z/fa)*(sfc-2*smg/(db*sin(sae)))'                 %球体垂直方向的受力平衡方程 对应公式2.16  
            
            F11=eval(F1); %对符号表达式求值；   
            F12=eval(F2);                  
            F21=eval(diff(F1,'sae'));     
            F22=eval(diff(F1,'sai'));         
            F31=eval(diff(F2,'sae'));                  
            F32=eval(diff(F2,'sai'));                     
            x1=[F11 F12];
            x2=[F21 F22;F31 F32];         
            
            s=inv(x2)*x1';
            x0=x0-s   
            eps=norm(s);
          end
         
    ae=x0(1)
    ai=x0(2)
    %%%布鲁和哈姆罗克借助最小二乘法用线性回归得到了第一类椭圆积分、第二类椭圆积分、椭圆率的下列简化方程:
    Rxe=(db*(dm+db*cos(ae)))/(2*dm);%%%%计算椭圆积分需要的中间变量
    Rxi=(db*(dm-db*cos(ai)))/(2*dm);
    Rye=(ffe*db)/(2*ffe-1);%%%%计算椭圆积分需要的中间变量
    Ryi=(ffi*db)/(2*ffi-1);
    eli1e=1.5277+0.6023*log(Rye/Rxe);  %第一类椭圆积分   
    eli1i=1.5277+0.6023*log(Ryi/Rxi);     
    eli2e=1.0003+0.5968*(Rxe/Rye);    %第二类椭圆积分   
    eli2i=1.0003+0.5968*(Rxi/Ryi);   
    elie=1.0339*(Rye/Rxe)^0.636;     %椭圆率参数，为接触椭圆长半轴与短半轴之比  
    elii=1.0339*(Ryi/Rxi)^0.636;   
    b=atan(sin(ae)/(cos(ae)+r));        %滚动体姿态角（自转轴线与x y坐标轴夹角）
    wb=-w/(((cos(ae)+tan(b)*sin(ae))/(1+r*cos(ae))+(cos(ai)+tan(b)*sin(ai))/(1-r*cos(ai)))*r*cos(b));%滚动体的自转速度   
    wm=w*(1-r*cos(ai))/(1+cos(ai-ae));%滚动体的公转速度   
    fc=(pi/12)*den*dm*db^3*wm^2*w;    %滚动体离心力
    mg=I*sin(b)*wb*wm;     %考虑滚动体相对于内套圈的陀螺力矩  
    sre1=db*cos(ae)/dm;   
    sri1=db*cos(ai)/dm;   
    Pe=(4-(1/ffe)+(2*sre1)/(1-sre1))/db;      %%球和外圈主曲率和     
    Pi=(4-(1/ffi)+(2*sri1)/(1-sri1))/db;      %%球和内圈主曲率和  
      
   
    ke=eli1e*((9/2/eli2e/Pe)*(1/pi/elie)^2)^(1/3); %滚动体和外套圈的接触刚度系数 对应公式2.24
    ki=eli1i*((9/2/eli2i/Pi)*(1/pi/elii)^2)^(1/3);   %滚动体和内套圈的接触刚度系数  对应公式2.24
    Qe=(fa/z+2*mg*cos(ae)/db)/sin(ae) %滚动体和外圈的接触力
    Qi=fa/(z*sin(ai))%滚动体和内圈的接触力
      
    deltai=Qi^(2/3)*ki;     %滚动体和内圈接触的变形与载荷的关系，对应公式2.28*P13
    deltae=Qe^(2/3)*ke;     %滚动体和内圈接触的变形与载荷的关系，对应公式2.29*P13
    deta=(re-0.5*db+deltae)*cos(sae)+(ri-0.5*db+ deltai)*cos(sai)-(ffi+ffe-1)*db*sin(a0);  

   
    Ke=1.5/eli1e*((9/2/eli2e/Pe)*(1/pi/elie/DE)^2)^(-1/3);% 和外圈接触刚度
    Ki=1.5/eli1i*((9/2/eli2i/Pi)*(1/pi/elii/DE)^2)^(-1/3);% 和内圈接触刚度
   
    KR=0;
    KA=0;
    KAN=0;
    for t=1:z   
        if isreal(Ke)==0 break;end
        if isreal(Ki)==0 break;end
        kri=Ki*cos(ai)^2;
        kai=Ki*sin(ai)^2;
        kre=Ke*cos(ae)^2;
        kae=Ke*sin(ae)^2;
        KR=KR+kri*kre*cos(2*pi*(t-1)/z)^2/(kri+kre);
        KA=KA+kai*kae/(kai+kae);   
        KAN=KAN+kai*kae*cos(2*pi*(t-1)/z)^2/(kai+kae);        
    end
KA;
KR;
KAN=KAN*dm*Ri;
Z=[KR KA KAN];

ChaChing

新手，不懂什么叫精华帖，哪里可以看到这些

ChaChing

这样是不是清楚一些呢

或许个人不善表诉, 致楼主可能还是没明白!?
楼主这程序又长且针对性专业, 楼主你认为有几人会细看? 若是你, 会吗?
11F不是已经点出楼主这程序第一个复数出现的位置了!
建议从此回头逐一检查吧
good luck

happy

ME! 发表于 2012-11-19 14:15

                               
登录/注册后可看大图

新手，不懂什么叫精华帖，哪里可以看到这些，
这样是不是清楚一些呢
function Z=dingyayujin7003(n,f)

参考http://forum.chinavib.com/forum. ... 94&fromuid=1691
看看是不是同类问题