2013-05-27 46 views
1

我有一个小问题。我有2个运动方程'ph'和'ph2'我不知道如何在x(1)> 0.111时设置ODE来停止计算'ph',然后开始再次计算'ph2'到0.111,之后在一张图上绘制'ph'+'ph2'取决于时间'w',我认为我必须设置一些时间限制,但不知道如何去做。我使用帮助但对我没有好处。MatLab ODE启动/停止条件

[t,y] = ode45(@ph,[0,w_max],[0,0]); 

function dx = ph(tt,x) 
global F1 c m_c Ff p w s ln f_t sig dstr Ren pn Fex Fzmax xz xn l Fz mn 
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-mn; 
end 
dx=[x(2);((spline(w,p,tt)*s)-Fpp-Ff-Fex-Fz)./m_c]; 


[t2,y2] = ode45(@ph2,[0,w_max],[0,0]); 

function dx=ph2(tt,x) 

    global Fv m_c g f alfa Fzp c m_nbp 

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

    if x(1)>0.44 

    m_c=m_c+m_nbp 

    end 

    dx = [x(2);((x(1)*c)-Ft-Fv)/m_c]; 
+0

你的意思是说你真的想建立一个分段函数作为ODE吗? – wakjah

+0

我不知道更好的解决方案:o/ – user2401142

回答

4

要停止计算pH值时x(1) > 0.111可以使用活动地点财产(关于如何使用它manual pageexample)。在实践中,它是在每个时间步骤评估的函数,如果返回的值为0,则ode45会停止积分。

添加功能

function [value,isterminal,direction] = events(t,y) 
% Locate the time when y passes through 0.111 in all 
% directions and stop integration. 
value = y(1) - 0.111; % Detect y=0.111 
isterminal = 1;  % Stop the integration 
direction = 0;   % All direction 

,并呼吁

[t,y] = ode45(@ph,[0,w_max],[0,0],options); 

之前options = odeset('Events',@events)设置但考虑到pH值的输出dx=[x(2); ...],为了做X上的一个检查(1)你需要输出也是这个变量 - 就像dx=[x(1); x(2);...]

希望这会有所帮助。

+0

这可能有助于查看Matlab附带的'ballode'演示。在命令窗口中输入'edit ballode'来查看代码。请注意,在这个例子中,他们如何使用前一个信息进行后续迭代,这可以大大提高效率:'options = odeset(options,'InitialStep',t(nt)-t(nt-refine),'MaxStep',t (NT)-t(1));'。 – horchler

+0

另外,使用全局变量不是一个好主意。您应该使用[匿名函数]传递参数(http://www.mathworks.com/help/matlab/matlab_prog/anonymous-functions.html)。查看我对[这个问题]的答案(http://stackoverflow.com/questions/16598558/finding-the-maxima-of-a-function-using-ode45/16639006#16639006),了解如何将参数传递给您的示例集成函数和使用'ode45'的事件函数。 – horchler

+0

Thx很多家伙。做得很好! 我需要工作,之后,我可以调试和优化。 – user2401142