func_1.m

Matlab数值解算法实现代码

%对数据进行后处理的函数，此函数中输出经过线形插值后的主凸轮和回凸轮的向径和压力角的值。
function func
[sol,kk]=fun;
sol=sol;                                       %输出没有对回凸轮的数据按凸进行处理以前的原始数据。
kk=kk;
%**************************************************************************
%**************************************************************************
%判定因子sequence用来决定是否对回凸轮的数据按凸轮转角fi进行排序或者按凸轮的向径角thita进行处理，
%sequence=1:按凸轮的向径角进行计算,sequence=2:按凸轮转角fi进行排序。
sequence=2;          
if sequence==1
    sol(360/kk+1,7)=sol(360/kk+1,7)+360;       %按凸轮转角fi对回凸轮的数据进行重新排序。
    for t=1:1:360/kk+1
        for s=t+1:1:360/kk+1
            if sol(t,7)>sol(s,7)
                m=sol(t,7:11);
                sol(t,7:11)=sol(s,7:11);
                sol(s,7:11)=m;
            end
        end
    end
else                                          %按凸轮的向径角进行计算。
    for w=2:1:360/kk+1
        if sol(w,7)<=sol(1,7)
            sol(w,7)=sol(w,7)+360;
        end
    end
    media=sol(1,7);
    for z=1:1:360/kk+1
        sol(z,7)=sol(z,7)-media;
    end
end        
sol=sol;                                     %输出处理后的结果
%对数据进行线形插值，得到所需的主凸轮和回凸轮向径和压力角的值。
a=0:kk:360;                                  %对主凸轮线性插值。
for dita=0:kk:360                            
    m=round(dita/kk+1);
    for v=1:1:360/kk
        while and(dita>=sol(v,2),dita<sol(v+1,2))
            jj1=sol(v,5)+((dita-sol(v,2))/(sol(v+1,2)-sol(v,2)))*(sol(v+1,5)-sol(v,5));
            kk1=sol(v,6)+((dita-sol(v,2))/(sol(v+1,2)-sol(v,2)))*(sol(v+1,6)-sol(v,6));
            break
        end
    end                                      %输出主凸轮插值后的结果。
    jj1=jj1;
    re(m,1)=dita;
    re(m,2)=jj1;
    re(m,3)=kk1;
end
re(1,4)=sol(1,10);                           %由于回凸轮数据的第一项不为零，无法迭代计算，所以首先对第一项进行计算。
re(1,5)=sol(1,11); 
for dita=kk:kk:360                           %对回凸轮线性插值。
    m=round(dita/kk+1);
    for v=1:1:360/kk
        while and(dita>=sol(v,7),dita<sol(v+1,7))
            jj2=sol(v,10)+((dita-sol(v,7))/(sol(v+1,7)-sol(v,7)))*(sol(v+1,10)-sol(v,10));
            kk2=sol(v,11)+((dita-sol(v,7))/(sol(v+1,7)-sol(v,7)))*(sol(v+1,11)-sol(v,11));
            break
        end
    end
    re(m,4)=jj2;                            %输出主凸轮插值后的结果。
    re(m,5)=kk2;
end
%**************************************************************************
%**************************************************************************
%输出主凸轮和回凸轮插值后的最后结果。
re=re
figure;                                    %输出主凸轮和回凸轮的轮廓线。
plot(sol(:,3),sol(:,4),sol(:,8),sol(:,9))
xlabel('X轴,mm');
ylabel('Y轴,mm');
hold on
figure                                     %输出主凸轮和回凸轮的向径值。
plot(re(:,1),re(:,2),re(:,1),re(:,4));
xlabel('向径角θ,Degree');
ylabel('向径ρ,mm')
hold on;
figure;                                    %输出主凸轮和回凸轮的压力角值。
plot(re(:,1),re(:,3),re(:,1),re(:,5));
xlabel('向径角θ,Degree');
ylabel('压力角α,mm');




%子程序fun是用来求出主凸轮和回凸轮的坐标值和压力角值的。子程序的接口为sol和kk，sol为%(360/ kk+1)*11的矩阵。kk为凸轮转角的间隔值。
function [sol,kk]=fun
%求出凸轮转角为0时D点的初始参数DY0和摆杆的初始角度a0。
format short;
kk=0.4;
l1=216;
l2=280;
l3=185;
l4=56;
bx=0;
by=0;
bx1=0;
by1=0;
dx=175;
a0=acos((l1^2+l3^2-l4^2)/(2*l1*l3));              %摆杆的初始转角a0。
cx0=l1*cos(a0);
cy0=l1*sin(a0);
a1=acos((cx0-dx)/l2);
dy0=cy0+l2*sin(a1);                           %D点的初始坐标DY0。
result=zeros(360/kk+1,7);
xx=zeros(360/kk+1,2);
yy=xx;
% *********************************************************************************
% *********************************************************************************
%循环开始求出凸轮的坐标值。
for fi=0:kk:360
    if and(fi>=0,fi<120)               %求出D点的位移和速度。
        s=25*(10*(fi/120)^3-15*(fi/120)^4+6*(fi/120)^5);
        s1=(25*3/(2*pi))*(30*(fi/120)^2-60*(fi/120)^3+30*(fi/120)^4);
    elseif and(fi>=120,fi<250)
        s=25;
        s1=0;
    elseif and(fi>=250,fi<330)
        s=25-25*(10*((fi-250)/80)^3-15*((fi-250)/80)^4+6*((fi-250)/80)^5);
        s1=-1*((25*9)/(4*pi))*(30*((fi-250)/80)^2-60*((fi-250)/80)^3+30*((fi-250)/80)^4);
    else
        s=0;
        s1=0;
    end
    [bita,bita1]=link(s,s1,dy0,a0);       %调用link子程序求出摆杆角度bita和角加速度bita1。
    res=cam(fi,bita,bita1);            %调用cam子程序求出凸轮在某一凸轮转角fi坐标值。
    result(fi/kk+1,:)=res;              %通过循环得到0到360度主凸轮和回凸轮的坐标值。
end
xx(:,1)=result(:,2);
xx(:,2)=result(:,5);
yy(:,1)=result(:,3);
yy(:,2)=result(:,6);
tht=zeros(360/kk+1,2);
for u=1:1:2                          %求出主凸轮和回凸轮的向径角θ。
    for t=1:1:360/kk+1
        if xx(t,u)>=0 
            if yy(t,u)>=0
                tht(t,u)=(atan(yy(t,u)/xx(t,u)))*180/pi;
            else
                tht(t,u)=(2*pi+atan(yy(t,u)/xx(t,u)))*180/pi;
            end
        else
            tht(t,u)=(pi+atan(yy(t,u)/xx(t,u)))*180/pi;
        end
    end
end
sol=zeros(360/kk+1,11); 
sol(:,1)=result(:,1);                              %输出主凸轮和回凸轮的转角φ。
sol(:,2)=tht(:,1);                                 %输出主凸轮的向径角θ。
sol(:,3)=result(:,2);                              %输出主凸轮的X坐标。
sol(:,4)=result(:,3);                              %输出主凸轮的Y坐标。
sol(:,5)=sqrt(result(:,2).^2+result(:,3).^2);      %输出主凸轮的向径值ρ。
sol(:,6)=result(:,4);                              %输出主凸轮的压力角α。
sol(:,7)=tht(:,2);                                 %输出回凸轮的向径角θ。
sol(:,8)=result(:,5);                              %输出回凸轮的X坐标。
sol(:,9)=result(:,6);                              %输出回凸轮的Y坐标。
sol(:,10)=sqrt(result(:,5).^2+result(:,6).^2);     %输出回凸轮的向径值ρ。
sol(:,11)=result(:,7);                             %输出回凸轮的压力角α。
sol=sol;
% **********************************************************************************
% **********************************************************************************
%LINK子程序求出摆杆的角位移dita和角加速度dita1。
function [bita,bita1]=link(s,s1,dy0,a0)
    bx=0;
    by=0;
    bx1=0;
    by1=0;
    dx=175;
    dx1=0;
    dy=dy0+s;
    dy1=s1;
    l1=216;
    l2=280;
    d=sqrt((dy-by)^2+(dx-bx)^2);
    dita=atan((dy-by)/(dx-bx));
    gama=acos((d^2+l1^2-l2^2)/(2*d*l1));
    bit=dita-gama;
    cx=l1*cos(bit);
    cy=l1*sin(bit);
    bita=bit-a0;
    bita1=((dx1-bx1)*(cx-dx)+(dy1-by1)*(cy-dy))/((cy-dy)*(cx-bx)-(cy-by)*(cx-dx));
% **********************************************************************************
% ********************************************************************************
% 计算凸轮的坐标值的子程序
function res=cam(fi,bita,bita1)
    fi=(fi/180)*pi;                   
    b=bita;                          
    b1=bita1;
    r01=130;
    r02=75.5;
    r1=21;
    r2=17.5;
    ll1=185;
    ll2=193;
    xy=zeros(1,5);
    res=zeros(1,5);
    for ii=1:1:2
        if ii==1
            l=ll1;
            r=r1;
            r0=r01;
            m=1;
            n=-1;
            k=1;
        else
            l=ll2;
            r=r2;
            r0=r02;
            m=-1;
            n=-1;
            k=1;
            ii=ii+2;
        end
        a=sqrt(130^2+185^2);
        bita=acos((a^2+l^2-r0^2)/(2*a*l));
        labu=atan((l-a*cos(bita))/(2*a*sin(bita)));
        af=atan((l*(1-n*k*b1)-a*cos(bita+k*b))/(a*sin(bita+k*b)));
        x=l*sin(labu+k*b-n*fi)-m*r*cos(labu+k*b-n*fi-af)-a*sin(labu-bita-n*fi);
        y=(l*cos(labu+k*b-n*fi)+m*r*sin(labu+k*b-n*fi-af)-a*cos(labu-bita-n*fi))*n;
        aa=abs(af*180/pi); 
        res(1,1)=fi/pi*180;
        res(1,ii+1)=x;
        res(1,ii+2)=y;
        res(1,ii+3)=aa;
    end