-1

我有一个 ODE 求解器运行良好且流畅，但我需要将所有内容绘制在一个图中。连接图（1）+（3）和图（2）+（4），我必须设置开始和停止条件，但这对我不起作用，我在死胡同。我正在尝试通过 x_m 设置结束条件，但没有结果。

options = odeset('Events',@events);

[t,y] = ode45(@ph1,[0,w_max],[0,0], options);
figure(1),plot(t,y(:,1));

x_n = y(:,1);
v_n = y(:,2);
x_m = x_n(end);
v_m = v_n(end);
q = max (t);


 d_v1 =diff(y(:,2));    
%d_t1 = diff(t);
%a_c1 = d_v1./d_t1;
 t_c1 = t(1:end-1);
%t_h1 = d_t1./2;

figure (2)
plot((t_c1),d_v1,'r') 
set(gca,'FontName','Times New Roman CE')
title('Rychlost')
xlabel('\it t\rm [s]')
ylabel('\it v_n\rm [m*s^{-1}]')
hold on

[t,y] = ode45(@ph2,[0,w_max],[0,0], options);
figure(3),plot(t,(y(:,1)));

d_v =diff(y(:,2));     
%d_t = diff(t);
%a_c = d_v./d_t;
t_c = t(1:end-1);
%t_h = d_t./2;

figure(4),plot((t_c),(d_v), 'g' );

% d_v2 =diff(d_v);     
% d_t2 = diff(d_t);
% a_c2 = d_v2./(d_t2.*d_t2);
% t_c2 = t(1:end-2);

% figure(5),plot((t_c2),a_c2 , 'r');

function [value,isterminal,direction] = events(t,y)

global ch

value = y(1) - ch;  
isterminal = 1;        
direction = 0;        




function dx = ph1(tt,x)
global F1 c m_c Ff p w s ln f_t sig dstr Ren pn Fex Fzmax xz xn l Fz m_n

Fpp = F1 + c*x(1); 

 if pn<0
     pn=abs(pn);
 end

if x(1)<ln

    pn=spline(w,p,tt)-((2*sig)/dstr*Ren);    
    Fex=3.1416.*f_t.*pn.*(ln-x(1));
end

if x(1)<42e-5
     Fz = Fzmax*(1-(1/xz)*(x(1)+l));   
end

if x(1)>44e-3
    m_c=m_c-m_n;
end


dx=[x(2);((spline(w,p,tt)*s)-Fpp-Ff-Fex-Fz)./m_c];

 function dx=ph2(tt,x)

 global Fv Ft m_z g f Fzp alfa m_nbp c

        Ft=m_z*g*f;
        Fv = 2*f*(Fzp/cos(alfa));

        if x(1)>0.44

        m_z=m_z+m_nbp

        end

        dx = [x(2);((x(1)*c)-Ft-Fv)/m_z];
4

2 回答 2

0

我用这个：

yy=y(:,1);
tt=t;
figure(3),plot(t+tt(end),abs((y(end:-1:1,1))),tt,yy);

其中 (t,y) 来自第一个 ODE， (tt,yy) 来自第二个。也许可以帮助某人

于 2013-05-30T06:35:13.387 回答
0

好的，让我先说一下：我不确定我是否有答案，因为我无法复制您的问题（您的代码无法运行......）并且由于某种原因我仍然无法发表评论。但我想我可以帮助你。

I. 似乎您正在分别求解两个 ODE，但两次都使用相同的时间 vecotr [0 w_max]。您可以（我建议）将两个 ODE 放入一个函数 dx = diffeq(args) 并简单地将 dx 设为行向量。您可能遇到的问题是，任何 ODE 求解器都会对不同的数学函数使用不同的时间步长。现在你可能想用“MaxStep”和“MinStep”来调节那个时间步长，但不要那样做。与简单的积分器相比，您将失去速度和准确性以及 matlab 中 ODE 求解器的所有优势。你应该把微分方程放到同一个函数中，即使它们在逻辑上根本没有联系。您可以更轻松地处理数据。

二、一旦您在同一个函数中求解两个微分方程，因此在 ode45/ode23 的同一运行中，它们将具有完全相同的时间向量。此外，解向量将具有相同的大小。那么比较它们应该是一件容易的事。但是，我不明白您如何无法“连接”这些数字。因为如果你总是用 'plot(x,y)' 而不是 'plot(y)' 绘图，那么它们应该完美对齐。例如：

 figure

 plot(t1,y1)

 hold on

 plot(t2, y2)

你会注意到我使用了“t1”和“t2”，你也应该这样做。在您当前的代码中，您正在用第二个“t”覆盖第一个“t”，这对调试不利。

三、只是另一个建议：你正在做类似的事情

 x_n = y(:,1);

 v_n = y(:,2);

相反，您可以为索引创建一个结构并将变量名称存储在其中。然后，您永远不会按数字访问变量，而只能通过这些索引访问变量。例如：

 ind.x_n = 1;

 ind.v_n = 2;

稍后在微分方程中： dy(ind.x_n) = ... dy(ind.v_n) = ... 一旦您拥有更多代码和更多状态，它将对您有很大帮助。

四。这让我想到了最后一点。更多州！你用全局变量做的事情非常非常糟糕！调试几乎是不可能的，并且您的代码会变慢。我在 diffeq 中使用了一次 ONE 全局变量，速度下降了 30% 左右。

  • 如果您有变量，您需要“生存”一个时间步，假设您想在下一步中使用旧的“F1”，那么您应该从中创建一个状态。这当然意味着您实际上可以推导出它。我不知道你在模拟什么，所以我不能说。去尝试一下！

  • 如果你有变量，你只在一个时间步计算，并且只在那个特定的时间步需要它们，那么真的没有必要使变量成为全局变量。只会害人！

我希望这不是太多，它对你有帮助。如果没有，请提供指向 git 存储库的链接，例如带有可执行代码的链接。这样就成了猜谜游戏了……

于 2013-05-28T18:18:19.693 回答