用ode45求解微分方程

方程是这样的y=dsolve('D2y-k*(r0-y)-m*g+c*Dy=0','D2y=9.81,y(0)=r0,Dy(0)=v0')但是求不出解析解,于是用ode4... 方程是这样的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的增函数,明显和阻尼运动不符。高手看看问题出在哪里呢?
展开
 我来答
liuliangsxd
2013-08-24 · TA获得超过1423个赞
知道小有建树答主
回答量:1174
采纳率:100%
帮助的人:651万
展开全部

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');

 

推荐律师服务: 若未解决您的问题,请您详细描述您的问题,通过百度律临进行免费专业咨询

为你推荐:

下载百度知道APP,抢鲜体验
使用百度知道APP,立即抢鲜体验。你的手机镜头里或许有别人想知道的答案。
扫描二维码下载
×

类别

我们会通过消息、邮箱等方式尽快将举报结果通知您。

说明

0/200

提交
取消

辅 助

模 式