常微分方程求解问题

boshi

有一个项目，模拟球自由落体与地面发生碰撞的过程。方程很简单，就是一个二阶常微分方程。难点在于要时刻判断什么时候与地面接触，接触之后经弹性力作用发生反弹，之后再发生碰撞直到静止。需要matlab编程，且用变步长数值解法。请各位帮忙.

chenzhian

帅哥，这样一个问题貌似可以直接可以求出方程的，没必要用数值分析的方法搞的哟。

xyz999

1# boshi

matlab的帮助文件里有你所说的例子。试试在matlab的command window输入“ballode”。

下面是matlab例子的源程序：
function ballode
%BALLODE  Run a demo of a bouncing ball.
tstart = 0;
tfinal = 30;
y0 = [0; 20];
refine = 4;
options = odeset('Events',@events,'OutputFcn', @odeplot,...
                 'OutputSel',1,'Refine',refine);
set(gca,'xlim',[0 30],'ylim',[0 25]);
box on
hold on;
tout = tstart;
yout = y0.';
teout = [];
yeout = [];
ieout = [];
for i = 1:10
  % Solve until the first terminal event.
  [t,y,te,ye,ie] = ode23(@f,[tstart tfinal],y0,options);
  if ~ishold
    hold on
  end
  % Accumulate output.  
  nt = length(t);
  tout = [tout; t(2:nt)];
  yout = [yout; y(2:nt,:)];
  teout = [teout; te];    % Events at tstart are never reported.
  yeout = [yeout; ye];
  ieout = [ieout; ie];
  ud = get(gcf,'UserData');
  if ud.stop
    break;
  end
  
  % Set the new initial conditions, with .9 attenuation.
  y0(1) = 0;
  y0(2) = -.9*y(nt,2);
  % A good guess of a valid first time step is the length of
  % the last valid time step, so use it for faster computation.
  options = odeset(options,'InitialStep',t(nt)-t(nt-refine),...
                           'MaxStep',t(nt)-t(1));
  tstart = t(nt);
end
plot(teout,yeout(:,1),'ro')
xlabel('time');
ylabel('height');
title('Ball trajectory and the events');
hold off
odeplot([],[],'done');
% --------------------------------------------------------------
function dydt = f(t,y)
dydt = [y(2); -9.8];
% --------------------------------------------------------------
function [value,isterminal,direction] = events(t,y)
% Locate the time when height passes through zero in a
% decreasing direction and stop integration.
value = y(1);     % Detect height = 0
isterminal = 1;   % Stop the integration
direction = -1;   % Negative direction only

boshi

你好，我把问题做成了一个文件付在附件中，请你指导一下，这类问题该如何编程求解呀？谢谢

xyz999

4# boshi

    既然楼主提供“有偿求助”，看来楼主不缺钱，那我就不客气啦！各位看官觉得多少钱合适呢？10.0元合适不？
    我提的条件：楼主如果碰见了捡垃圾的且年龄大于60岁的老人的话，就把这10元钱捐给他（她）吧，如果生活不艰难，他们是不会出来捡垃圾的。

    言归正传，楼主的问题就是一个“事件判断”的问题，matlab已经提供了类似的功能，下面是源程序。程序中的m, n, k, ce, h, r等参数是我随便设置的，可根据实际情况做相应的更改。

function [tout yout teout] = myfun()
% calculate the trajectory of a falling ball
% tout --- time
% yout --- trajectory of the ball
% teout --- time when the ball touch the floor

m = 2.0; % mass of the ball
n = 2.0;
k = 2e5; % contact stiffness
ce = 0.6; % Coefficient of restitution
h = 5; % initial height of the ball
r = 0.04; % radius of the ball

stepn = 20; % maximum bouncing times
ystop = 1e-6; % if disp<ystop, and velocity < vstop, stop simulation
vstop = 1e-6;

tstart = 0;
dt = 10; % time for integration
y0 = [h; 0.0]; % initial value
optn = odeset('RelTol',1e-5,'AbsTol',1e-8,'Events',@eveq);

detvn = 1.0; % initial impact velocity
tout = tstart;
yout = y0';
teout = [];
yeout = [];
ieout = [];

% final interaction
yend = (m*9.8/k)^(1/n);

for ii = 1:stepn
[t,y,te,ye,ie] = ode45(@beq,[tstart tstart+dt],y0,optn);

% Accumulate output.
nt = length(t);
tout = [tout; t(2:nt)];
yout = [yout; y(2:nt,:)];
teout = [teout; te];
yeout = [yeout; ye];
ieout = [ieout; ie];

% judge whether to stop
if (abs(y(nt,1)-(r-yend))<ystop)&&(abs(y(nt,2))<vstop)
break;
end

if ~isempty(ie)
detvn = y(nt,2);
end

tstart = t(nt);
y0 = y(nt,:)';
end
plot(tout,yout(:,1),'b-',teout,yeout(:,1),'ro')
xlabel('time');
ylabel('height');
grid on;
title('Ball trajectory');

% motion equation
function dydt = beq(t,y)
delty = r - y(1);
if delty < 0
% no contact
dydt = [y(2); -9.8];
else
ff = k/m * delty^n * (1.0 + 0.75*(1-ce*ce) * y(2)/detvn) - 9.8;
dydt = [y(2); ff];
end
end

% event equation
function [value,isterminal,direction] = eveq(t,y)
value = (r - y(1));
isterminal = 1;
direction = 1;
end
end



程序输出参数：
tout：时间
yout：小球的轨迹
teout：小球与地面开始碰撞的时间

chenzhian

5# xyz999
精彩啊！xyz兄，评论精彩，解答也精彩！对于这个问题，一个平台，一个球，做个动画挺好的，也很直观。哈哈！建议哈。

boshi

谢谢您精彩的回复，但是还有些问题需要请教您，比如碰撞速度是需要判断的，而不是定值（detvn = 1.0; % initial impact velocity）
qq229612352,希望可以交流一下。

xyz999

7# boshi

detvn 不是定值，在程序中会改变的，我给的那个值只是用来初始化detvn这个变量而已。
if ~isempty(ie)
   detvn = y(nt,2); % 在这里设置的
end

boshi

clc;
clear all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%碰撞接触刚度计算
r2=0.02;      %球半径；
r3=0.04;      %假设地面也为一球面
v2=0.3;       %泊松比
v3=0.3;       %泊松比
E2=207000000000;%弹性模量
E3=207000000000;%弹性模量
sg2=(1-v2^2)/E2;
sg3=(1-v3^2)/E3;
k=4/(3*(sg2+sg3))*(r2*r3/(-r2+r3))^0.5;%接触刚度。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=1.5;
e=0.9;            %e恢复系数；
m=10;             %球质量kg
g=9.8;            %重力加速度
h=1;              %自由落体初始高度m
vvcc=[];
q=[];
FF=[];
VVn=[];
dd=[];
y0=[h;0];          %初值（位置/速度）
tt=0:0.00001:4;    %计算时间，注意步长选取很小，否则计算的反弹力会很大，以至球反弹的高度比自由落体前的高度还要高，不符合实际（问题所在）。
for j=1:400000     %（按时间划分为100000个计算周期）
    d=y0(1)+r2;    %球最下端位置
    Vn=y0(2);      %接触相对速度
    if d<=0
        dd=[dd,d];        %记录球最下端位置
        vvcc=[vvcc,Vn];   %同时，记录球的速度
        vcc=vvcc(1);      %取出碰撞接触开始时刻速度；
        F=k*(abs(d))^n*(1+3*(1-e^2)/4*Vn/vcc);  %计算接触反弹力，与前期计算的刚度K和接触深度abs(d)及速度有关。   
        FF=[FF,F];        %记录反弹力的变化
    else
        F=0;        %没有发生接触，所以接触力为零
        vvcc=[];    %清空速度记录
    end
    Fn=F;           %取出反弹力，用于以下微分方程求解
    t0=tt(j);tf=tt(j+1);%开始计算
    tspan=linspace(t0,tf);
    [t,y]=ode45(@fangcheng,tspan,y0,[],m,Fn,g);
    y0=y(size(y,1),:); %y的最后一行所有列付给y0。也就是将本次计算结果，作为下次计算的初值。
    q=[q;y0];%收集结果
end
plot(tt(1:400000),q(:,1));%绘制球自由落体运动


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Dydt=fangcheng(t,y,m,Fn,g);
Dydt=[y(2);
      (-m*g+Fn)/m];


上面是我编写的程序，请大家指导

%问题描述：为了保证接触反弹力计算结果的正确性，计算中步长选取很小，导致计算时间很长。
%程序修改方向：在接触之前采用大步长计算，在即将发生接触之前，步长马上调整为小步长。（问题所在）
计算结果见附件。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xyz999

如果不想使用matlab提供的event方法，楼主可以朝着下面的方向思考解决方法：

t0=tt(j);tf=tt(j+1);%开始计算
在这一行插入一个判断，判断y的大小和速度方向。如果：(下面用的的abc为一人为设定的临界值）
    1. y>abc and y' > 0，说明小球已经离开地面，并且向上走，可以使用最大的步长，dt = ...；
    2. y>abc and y' < 0，说明小球正在冲向地面，但还有一段距离，可以使用较大的步长，dt = ...；
    3. y<abc，在接触区附近，使用较小的步长，dt = ...。
tspan=[t0    t0+dt];
[t,y]=ode45(@fangcheng,tspan,y0,[],m,Fn,g);
y0=y(size(y,1),:); ........

没有实际做，仅仅是思路而已。

boshi

谢谢您的指导。