微分方程编程求解（有偿求助）
有一个项目，模拟球自由落体与地面发生碰撞的过程。方程很简单，就是一个二阶常微分方程。难点在于要时刻判断什么时候与地面接触，接触之后经弹性力作用发生反弹，之后再发生碰撞直到静止。需要matlab编程，且用变步长数值解法。请各位帮忙，有偿酬谢！该不是什么课的大作业吧……{:{01}:} 回复 2 # Rainyboy 的帖子
请问，你会求解这类问题吗？ 初步分析了下，这个问题并不复杂。没有碰撞的时候是自由落体运动，可以用大步长，出现碰撞时用小步长进行数值计算，判断条件题目已给。变步长的数值计算方法应该是比较成熟的了。自己找找看
分析问题的过程需要自己进行，中间有问题可以询问。整个问题放在这里，不加自己的分析和思考那是不太负责的 本帖最后由 greatchina 于 2010-11-24 17:44 编辑
回复 4 # hustxyon
你好，这是我已经编好的程序，可以解决问题，但是计算的步数太多了，希望可以修改一下步长。请多多指教，谢谢
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];
%问题描述：为了保证接触反弹力计算结果的正确性，计算中步长选取很小，导致计算时间很长。
%程序修改方向：在接触之前采用大步长计算，在即将发生接触之前，步长马上调整为小步长。（问题所在）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
直接用ode45求解器求解是实现不了变步长的，有两个思路可选择
1、常微分方程变步长数值算法已经有很多成熟的算法，可以解决变步长的问题，你可以去查找一下相关资料。
2、程序里用while循环其实就可以用ode45实现两个状态用两种步长的功能，再想想看 请问这个问题解决了吗？我也在做这方面的研究，急切需要高手指点，万分感谢！
页:
[1]