[已答复] ODE求解非线性二元二阶常微分方程(动力学微分方程)

[复制链接]
free北宸 发表于 2017-11-10 10:33:41
请大家帮我看个问题可好。
我想用ode解一个变系数二元二阶常微分方程。
这个常微分方程的系数是随时间变化的,分成了四种情况,我需要每一个步长检验一下系数矩阵是不是变化了。
我写了一个程序,自己定义一个大步长dt,该步长内用ode算,每dt查一下K是否变化,就是出来的结果不是简谐,与原始文献给出的结果对不上,复现不了别人的结果,我这个课题的思路也就推进不下去。
课题卡住两个星期了,请大家有时间的话可不可以帮我看看,谢谢!
下面是方程:









下面是我自己编写的程序:
m=1*10^3;
J=2.5*10^3;
b=2;%初始参数
dt=0.0001;%在ode外人为设定适用步长
Xt=zeros(10000,4);%将每一步长计算得到的解累加
Tt=zeros(10000,1);%将每一步长的ode步长累加
T_row=0;%时间T-total行数初始时为0
x0=[0,0,0,0.5];%初始条件u、u'、θ、θ'
for n=1:100000
    t=n*dt;
    x=x0;%先给x赋初值
if x(1)+(b/2)*x(3) > 0 || x(1)+(b/2)*x(3) == 0;
    k1=1*10^6;
else  k1=2*10^7;
end%判断拉压刚度不同弹簧1刚度值
if x(1)-(b/2)*x(3) > 0 || x(1)-(b/2)*x(3) == 0;
    k2=1*10^6;
else k2=2*10^7;
end%判断拉压刚度不同弹簧2刚度值   
[url=mailto:dx=@(t,x)[x(2]dx=@(t,x)[x(2[/url]); (-(k1+k2)*x(1)-(b/2)*(k1-k2)*x(3))/m; x(4); ((b/2)*(k1-k2)*x(1)+((b^2)/4)*(k1+k2)*x(3))/J];
%二元二阶常微分方程表达式,这里x1=u,x2=u',x3=θ,x4=θ'
[t,x]=ode45(dx,[(n-1)*dt,t],x0);%调用ode求解器解常微分方程
[row,column]=size(x);%[行,列]=x解得矩阵的维数
x0=x(row,:);%上一步长dt结束时的解是下一步长ode计算时的初始条件
for i=1:row-1
    Tt(T_row+i,1)=t(i,1);
    Xt(T_row+i,:)=x(i,:);
end%循环的目的是把人为定义步长dt内ode自己划分的步长累加到总的时间矩阵T-total中
T_row=T_row+row-1;
end
plot(Tt,Xt(:,4));%绘图u-t
% plot(Tt,Xt(:,3));%绘图θ-t

2 条回复


free北宸 发表于 2017-11-10 15:35:26
RS4QWO29$)H]6OOC]`QAP{3.png
抱歉刚才方程未能上传,这里补发。

无名居士 发表于 7 天前
朋友,您好。我也遇到和你相似的问题,不知道您最后如何解决了这个问题?麻烦可以分享一下吗?
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表