[求助] 求帮助！！关于利用matlab模拟二自由度含间隙碰撞系统的问题，本科毕业论文！！！！！

各位论坛的大神，小弟本科毕业论文，老师给的模型是这样一个含间隙的二自由度碰撞系统，研究它的分叉与混沌：左右各两个小球，各自连接着一个非线性的弹簧和阻尼器，加上激励力 具体的我上传了图片，还有附件。 问题描述完毕，我现在遇到的问题是：要在matlab里考虑这两个小球的碰撞，老师说要把时间步设置为1e-6pi，计算时间为500pi，并且去掉前面200pi的时间。但是！这个根本算不了哇，我这个电脑算这么多的话得算到猴年马月去了，然后我自就把时间缩短到50pi，时间步改为1e-3pi，但是计算起来也费时间，我现在就想问问各位： 1.这个问题必须要这么长的时间吗？有什么办法可以缩减一下？ 2.怎样判断这两个小球碰撞上了？ 我写的代码是：if abs(delta-abs(x(1,n+1)-x(2,n+1)))<=10^(-2)   delta是间隙的距离。x（1，n+1）是第n+1步时第一个球的位置，x（2，n+1）是第n+1步时第二个球的位置。不知道这样写对不对。 昨晚算的结果我截了图，传附件了。 希望大哥们能帮小弟看看。 对了，还有我自己编的matlab代码： 第一个dyfun函数，用于存放控制方程 function f=dyfun(x,y) f(1)=y(3); f(2)=y(4); f(3)=-0.01*y(3)*0.1-y(1)*1^2-0.01*y(1)^3+10*0.01*cos(x*2); f(4)=-0.01*y(4)*0.15-y(2)*1.5^2-0.012*y(2)^3+3*0.01*cos(x*2); f=f( : ); end -------------------------------------------------------------------------------------------------- 第二个函数，自己编写的4阶R-K方法，因为ode45里面我不知道怎么添加判断条件，判断两个小球是否撞上了。 function[t,x]=narunge4s(dyfun,tspan,x0,h) delta=input('delta='); t=tspan(1):h:tspan(2); x=zeros(length(x0),length(t)); x(:,1)=x0(: ); for n=1: (length(t)-1)     k1=feval(dyfun,t(n),x(:,n));     k2=feval(dyfun,t(n)+h/2,x(:,n)+h/2*k1);     k3=feval(dyfun,t(n)+h/2,x(:,n)+h/2*k2);     k4=feval(dyfun,t(n+1),x(:,n)+h*k3);     x(:,n+1)=x(:,n)+(h/6)*(k1+2*k2+2*k3+k4); %下面这三行是由碰撞过程中动量守恒和恢复系数得到的 %        lanbda=n1/n2 此处取做1 %        q+p=x(3,n+1)+x(4,n+1); %        q-p=-0.8*(x(3,n+1)-x(4,n+1));       if abs(delta-abs(x(1,n+1)-x(2,n+1)))<=10^(-2)    %这就是判断是否碰撞的语句，不知道对不对        q=(x(3,n+1)+x(4,n+1)-0.8*(x(3,n+1)-x(4,n+1)))/2;        p=(x(3,n+1)+x(4,n+1)+0.8*(x(3,n+1)+x(4,n+1)))/2;        x(3,n+1)=q;        x(4,n+1)=p;     这就是判断是否碰撞的语句，不知道对不对 %        stop=input('stop=');       end       n % 这是为了一边计算一边看图，隔100个点画一次    if mod(n,100)==0        figure(1);     plot(t(1:n),x(1,1:n),'r',t(1:n),x(2,1:n),'r-.');     legend('x(1)','x(2)');     hold on     xlabel('x(1)');     ylabel('x(2)');     figure(2);     plot(x(1,1:n),x(3,1:n),'r');     xlabel('x(1)');     ylabel('dx(1)/dt');     title('The 1st')     figure(3);     plot(x(2,1:n),x(4,1:n),'r-')     xlabel('x(2)');     ylabel('dx(2)/dt');     title('The 2nd')     hold on     end end ------------------------------------------------------------------------------------------------- 主程序，调用之前的计算。 clc clear all [t,x]=marunge4s(@dyfun,[0 50*pi],[0 0.5 0 0],10^-4*pi); 论文模型.PNG 两个小球的位移-时间图.PNG 第二个小球的位移-速度相图.PNG

回复此楼