我的程序主要是对碰撞过程弹簧阻尼系统的部分改写。这是一个求解微分方程并调用的程序，我想在主程序中调用子程序求解过程的结果，进而在主程序中得到相应的所求结果。但是我编的遇到的问题是主程序中调用的 Kl 显示等于零，跟实际不符。请大家给些修改意见，我自己一直改不好，很是难受。
求解过程程序：
function  dx=ImpactFn(t,x)
dx=zeros(2,1);
m=1;           
n=1.5;         
Ce=0.9;        
R1=10;        
R2=9.9;      
c=R1-R2;
E=2.07*10^5;  
v0=500;        
z=0.29;      
E1=E/(2*(1-z^2));
Ks=(pi*E1/8)*(2*x(1)*(3*c+2*x(1))^2/(c+x(1))^3)^(1/2);
p1=7700*10^(-9);
c1=5900*10^3;   
p2=863*10^(-9);  
c2=1530*10^3;   
Aa=125*0.25;     
L=50;         
Z1=p1*c1;
Z2=p2*c2;
B=25*(pi/180);   
B0=5*(pi/180);   
ht=0.1;           
e1=acos(c/(c+x(1)));   
l=2*R1*e1;      
A0=l*L;         
for i=1:length(e1)
    if e1(i)>B
         q(i)=B/B0;
    else q(i)=floor(e1(i)/B0);
    end

           Kl(i)=0         
        for j=1:q(i)
          hl(i)=c+(x(1)+c)*cos(pi-B0*j)+ht;
          Kl0(i)=Z1*c2*Aa/((Z1/Z2-Z2/Z1)*hl*(B/B0))
          Kl(i)=Kl(i)+Kl0(i);
        end
end
save shuju1 Kl
Kn=Ks*(((A0)-Aa*q/(B/B0))/(A0+10^(-15)))+Kl;      
dx(1)=x(2);                        
dx(2)=(-1/m)*Kn*x(1)^n*(1+(8*(1-Ce))*x(2)/(5*Ce*v0));
end
主程序：
clear;clc;
y0=[0;500];         
[t,x]=ode45(@ImpactFn,[0, 2.44e-4],y0);  %求解微分方程的调用格式  
n=1.5;        
Ce=0.9;                           
R1=10;         
R2=9.9;      
c=R1-R2;
E=2.07*10^5;  
v0=500;      
z=0.29;        
E1=E/(2*(1-z^2));
Ks=(pi*E1/8)*(2*x(:,1).*(3*(R1-R2)+2*x(:,1)).^2/(R1-R2+x(:,1)).^3)^(1/2);
[~,col,~]=find(Ks,1,'first');
Ks=Ks(:,col);
p1=7700*10^(-9);
c1=5900*10^3;   
p2=863*10^(-9);  
c2=1530*10^3;   
Aa=125*0.25;     
L=50;            
Z1=p1*c1;
Z2=p2*c2;
B=25*(pi/180);   
B0=5*(pi/180);   
ht=0.1;           
e1=acos(c./(c+x(:,1)));   
l=2*R1*e1;      
A0=l*L;      

for i=1:length(e1)
        if e1(i)>B
           q(i)=B/B0;     
          else q(i)=floor(e1(i)/B0);
        end   
end
load shuju1 Kl
Kl=Kl.';            
q=q.';
Kn=Ks.*(((A0+10^(-15))-Aa*q/(B/B0))./(A0+10^(-15)))+Kl;      
Fn=Kn.*x(:,1).^n.*(1+(8.*(1-Ce)).*x(:,2)/(5*Ce*v0));
figure(1);         
plot(t,x(:,1),'red');
xlabel('时间/s','Fontname', '宋体','FontSize',13);
ylabel('变形量/mm','Fontname', '宋体','FontSize',13);
hold on            
figure(2);
plot(t,x(:,2),'red');
xlabel('时间/s','Fontname', '宋体','FontSize',13);
ylabel('接触变形速度/mm·s-1','Fontname', '宋体','FontSize',13);
hold on                          
figure(3);           
plot(t,Fn,'red');   
xlabel('时间/s','Fontname', '宋体','FontSize',13);
ylabel('接触力/N','Fontname', '宋体','FontSize',13);
hold on            
figure(4);         
plot(x(:,1),Fn,'red');
xlabel('变形量/mm','Fontname', '宋体','FontSize',13);
ylabel('接触力/N','Fontname', '宋体','FontSize',13);
hold on            
figure(5);         
plot(x(:,1),Kn,'red');
xlabel('变形量/mm','Fontname', '宋体','FontSize',13);
ylabel('接触刚度/N*mm-1','Fontname', '宋体','FontSize',13);
hold on 返回小木虫查看更多