常微分方程求解问题
本帖最后由 boshi 于 2010-11-27 08:37 编辑有一个项目，模拟球自由落体与地面发生碰撞的过程。方程很简单，就是一个二阶常微分方程。难点在于要时刻判断什么时候与地面接触，接触之后经弹性力作用发生反弹，之后再发生碰撞直到静止。需要matlab编程，且用变步长数值解法。请各位帮忙. 帅哥，这样一个问题貌似可以直接可以求出方程的，没必要用数值分析的方法搞的哟。 1# boshi
matlab的帮助文件里有你所说的例子。试试在matlab的command window输入“ballode”。
下面是matlab例子的源程序：
function ballode
%BALLODERun a demo of a bouncing ball.
tstart = 0;
tfinal = 30;
y0 = ;
refine = 4;
options = odeset('Events',@events,'OutputFcn', @odeplot,...
'OutputSel',1,'Refine',refine);
set(gca,'xlim',,'ylim',);
box on
hold on;
tout = tstart;
yout = y0.';
teout = [];
yeout = [];
ieout = [];
for i = 1:10
% Solve until the first terminal event.
= ode23(@f,,y0,options);
if ~ishold
hold on
end
% Accumulate output.
nt = length(t);
tout = ;
yout = ;
teout = ; % Events at tstart are never reported.
yeout = ;
ieout = ;
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 = ;
% --------------------------------------------------------------
function = 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 你好，我把问题做成了一个文件付在附件中，请你指导一下，这类问题该如何编程求解呀？谢谢 4# boshi
既然楼主提供“有偿求助”，看来楼主不缺钱，那我就不客气啦！各位看官觉得多少钱合适呢？10.0元合适不？
我提的条件：楼主如果碰见了捡垃圾的且年龄大于60岁的老人的话，就把这10元钱捐给他（她）吧，如果生活不艰难，他们是不会出来捡垃圾的。
言归正传，楼主的问题就是一个“事件判断”的问题，matlab已经提供了类似的功能，下面是源程序。程序中的m, n, k, ce, h, r等参数是我随便设置的，可根据实际情况做相应的更改。
function = 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 = ; % 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
= ode45(@beq,,y0,optn);
% Accumulate output.
nt = length(t);
tout = ;
yout = ;
teout = ;
yeout = ;
ieout = ;
% 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 = ;
else
ff = k/m * delty^n * (1.0 + 0.75*(1-ce*ce) * y(2)/detvn) - 9.8;
dydt = ;
end
end
% event equation
function = eveq(t,y)
value = (r - y(1));
isterminal = 1;
direction = 1;
end
end
程序输出参数：
tout：时间
yout：小球的轨迹
teout：小球与地面开始碰撞的时间 5# xyz999
精彩啊！xyz兄，评论精彩，解答也精彩！对于这个问题，一个平台，一个球，做个动画挺好的，也很直观。哈哈！建议哈。 谢谢您精彩的回复，但是还有些问题需要请教您，比如碰撞速度是需要判断的，而不是定值（detvn = 1.0; % initial impact velocity）
qq229612352,希望可以交流一下。
7# boshi
detvn 不是定值，在程序中会改变的，我给的那个值只是用来初始化detvn这个变量而已。
if ~isempty(ie)
detvn = y(nt,2); % 在这里设置的
end
本帖最后由 boshi 于 2010-11-24 14:15 编辑
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=; %初值（位置/速度）
tt=0:0.00001:4; %计算时间，注意步长选取很小，否则计算的反弹力会很大，以至球反弹的高度比自由落体前的高度还要高，不符合实际（问题所在）。
for j=1:400000 %（按时间划分为100000个计算周期）
d=y0(1)+r2; %球最下端位置
Vn=y0(2); %接触相对速度
if d<=0
dd=; %记录球最下端位置
vvcc=; %同时，记录球的速度
vcc=vvcc(1); %取出碰撞接触开始时刻速度；
F=k*(abs(d))^n*(1+3*(1-e^2)/4*Vn/vcc);%计算接触反弹力，与前期计算的刚度K和接触深度abs(d)及速度有关。
FF=; %记录反弹力的变化
else
F=0; %没有发生接触，所以接触力为零
vvcc=[]; %清空速度记录
end
Fn=F; %取出反弹力，用于以下微分方程求解
t0=tt(j);tf=tt(j+1);%开始计算
tspan=linspace(t0,tf);
=ode45(@fangcheng,tspan,y0,[],m,Fn,g);
y0=y(size(y,1),:); %y的最后一行所有列付给y0。也就是将本次计算结果，作为下次计算的初值。
q=;%收集结果
end
plot(tt(1:400000),q(:,1));%绘制球自由落体运动
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Dydt=fangcheng(t,y,m,Fn,g);
Dydt=[y(2);
(-m*g+Fn)/m];
上面是我编写的程序，请大家指导
%问题描述：为了保证接触反弹力计算结果的正确性，计算中步长选取很小，导致计算时间很长。
%程序修改方向：在接触之前采用大步长计算，在即将发生接触之前，步长马上调整为小步长。（问题所在）
计算结果见附件。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 如果不想使用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=;
=ode45(@fangcheng,tspan,y0,[],m,Fn,g);
y0=y(size(y,1),:); ........
没有实际做，仅仅是思路而已。 谢谢您的指导。
页:
[1]