24小时热门版块排行榜

返回列表

装沉默

铁虫 (初入文坛)

应助: 0 (幼儿园)
金币: 13
帖子: 24
在线: 5.5小时
虫号: 3743473
注册: 2015-03-16
专业: 固体力学

[求助] 求帮助！！关于利用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