Question

用ode45求解微分方程

方程是这样的y=dsolve('D2y-k*(r0-y)-m*g+c*Dy=0','D2y=9.81,y(0)=r0,Dy(0)=v0')
但是求不出解析解，于是用ode45求解析解
定义的函数如下：
function dydt=meme(t,y) %y,y的一阶导，二阶导变量分别用y(1),y(2),y(3)表示
k=5e4*0.5
r0=0.0125
Density=2650

m=3.1416*R^3*Density*4/3
g=9.81
cc=2*sqrt(k*m)
b=0.5
c=cc*b
h0=1
v0=sqrt(2*(h0-R)*g);
dydt=[y(3);y(2);(m*g+k*r0-y(3)-c*y(2))/k];
程序如下：
t_end=.9;
x0=[g;v0;r0];
[t,y]=ode45('meme',[0,t_end],x0);
plot(t,y(:,2));
xlabel('t'); ylabel('y');
但是解出来的y(3),y(2)都是t的增函数，明显和阻尼运动不符。高手看看问题出在哪里呢？

Answer 1

clear all

clc

k=5e4*0.5;
r0=0.0125;
Density=2650;

m=3.1416*r0^3*Density*4/3;
g=9.81;
cc=2*sqrt(k*m);
b=0.5;
c=cc*b;
h0=1;
v0=sqrt(2*(h0-r0)*g);
dydt=@(t,y)([y(2);m*g+k*(r0-y(1))-c*y(2)]);
t_end=.9;
x0=[g;v0];
[t,y]=ode45(dydt,[0,t_end],x0);

d2y=m*g+k*(r0-y(:,1))-c*y(:,2);
subplot(1,3,1),plot(t,y(:,1)),xlabel('t'),legend('y');
subplot(1,3,2),plot(t,y(:,2)),xlabel('t'),legend('dy/dt');
subplot(1,3,3),plot(t,d2y),xlabel('t'),legend('d2y/dt');