[已答复] ode45解动力学微分方程,if语句怎么判断

[复制链接]
Lab小白 发表于 2021-4-8 10:36:00
本帖最后由 Lab小白 于 2021-4-8 11:31 编辑

我在用ode45解微分方程时,需要在function中加一个if语句控制变量,但却总是起不到控制效果,图像结果是以0.2计算的,哪个大哥可以帮忙看看,可以包网费:):)  

%%%%%% 这是我的function文件
function dy=myfun(t,y,~,w,omega,gamma,A,B,alpha,beta,m)
x=y(1);          % x-位移
v=y(2);          % v-速度
z=y(3);          % z-恢复力
dy(1,1)=v;
    if   v*(v+B*w*cos(w*t))>0  
dy(2,1)=-2*0.2*v-gamma*x*(omega^2)-(1-gamma)*z/m+B*(w^2)*sin(w*t);     
   else
dy(2,1)=-2*0.1*v-gamma*x*(omega^2)-(1-gamma)*z/m+B*(w^2)*sin(w*t);
   end
dy(3,1)=A*v-alpha*abs(v)*z-beta*v*abs(z);
end

%%%%%%这是命定文件
h_opt=odeset;
omega=1;         %参数大小
gamma=0.5;
A=1;
B=0.01;
alpha=0.6;
beta=0.2;
m=200;
w=0.75;

y0=[0,0,0];      %初值
tspan=0:0.01:200;
[t,y]=ode45('myfun',tspan,y0,h_opt,w,omega,gamma,A,B,alpha,beta,m);
plot(t,y(:,1))



1 条回复


TouAkira 发表于 2021-4-8 12:21:41
其实没有问题,所谓的"总是起不到控制效果"是错觉。

原因是 B = 0.01; w = 0.75; 这两个数值太小了,导致你的判断条件 v*(v+B*w*cos(w*t)) 是否为正,起作用的区间非常非常地狭窄。具体地讲,B, w, t全是常数,余弦函数本身值域在[ -1, 1 ]以内,于是该判断条件的两个零点分别是 v = 0,以及 v = -B * w * cos( w * t ),大致可以估计出只有当 v 在[ -B * w, B * w ]也就是[ -0.075, 0.075 ]这个区间以内(实际上由于余弦部分导致这个零点区间应该会更小)时,才会走else对应的分支。

另外就是你的两个分支,系数差别也不是很大,有两个控制分支的微分方程组,与没有分支的方程组,数值解的差别非常小,再加上你的绘图时间比较长,这些许偏差肉眼看不到。
  1. B = 0.01; w = 0.75;
  2. Para_1 = 0.5; Para_2 = 0.01; % 两个控制分支的系数,放在这里改起来方便,你的原始值是Para_1 = 0.2; Para_2 = 0.1;
  3. % 你的有两个控制分支的微分方程组
  4. NumericalSolution47 = @( t, y ) [ y( 2 );
  5.     ( y( 2 ) * ( y( 2 ) + B * w * cos( w * t ) ) > 0 ) .*  ( -2 * Para_1 * y( 2 ) - gamma * y( 1 ) * omega^2 - ( 1 - gamma ) * y( 3 ) / m + B * w^2 * sin( w * t ) ) + ...
  6.     ( y( 2 ) * ( y( 2 ) + B * w * cos( w * t ) ) <=  0 ) .* ( -2 * Para_2 * y( 2 ) - gamma * y( 1 ) * omega^2 - ( 1 - gamma ) * y( 3 ) / m + B * w^2 * sin( w * t ) );
  7.     A * y( 2 ) - alpha * abs( y( 2 ) ) * y( 3 ) - beta * y( 2 ) * abs( y( 3 ) ) ];
  8. % 只有上面控制分支的微分方程组,作为对比
  9. NumericalSolutionNoBranch = @( t, y ) [ y( 2 );
  10.     -2 * Para_1 * y( 2 ) - gamma * y( 1 ) * omega^2 - ( 1 - gamma ) * y( 3 ) / m + B * w^2 * sin( w * t ) ;
  11.     A * y( 2 ) - alpha * abs( y( 2 ) ) * y( 3 ) - beta * y( 2 ) * abs( y( 3 ) ) ];
  12. y0 = [0,0,0];      %初值
  13. tspan = 0 : 1e-2 : 40; % 把绘图时间跨度缩短,这样才明显
  14. [ t, y ] = ode45( NumericalSolution47, tspan, y0, h_opt );% 你的有两个控制分支的微分方程组
  15. [ t_n, y_n ] = ode45( NumericalSolutionNoBranch, tspan, y0, h_opt );% 只有上面控制分支的微分方程组,作为对比
  16. plot( t, y( :, 2 ), 'k-.', t_n, y_n( :, 2 ), 'r-', 'LineWidth', 0.5 );
复制代码

这样再看v的图像,两者之间的差别就很明显了,图我懒得发了,你自己运行对比就是了
您需要登录后才可以回帖 登录 | 注册

本版积分规则

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