Need help to stop the ode when v(2) is zero using ode events. I have tried in so many ways but couldn't get result.

1 次查看（过去 30 天）
Dereje
Dereje 2017-12-30
评论： Dereje 2018-1-3
%%Data for desnity with respect to depth
z = [2 3 5 7 10 15 20 25 30 40 50 60 70 80 90 100 125 150 160 175 200 225 250 275 300 325 350 375 400];
rho = [17.2731684 17.1649375 21.43455647 22.4140625 23.86332207 24.3746967 24.70487685 24.6003125 24.8933125 25.42772826 26.03220776 26.439625 26.8151875 26.86830797 27.1949375 27.34406944 27.5551875 27.728625 27.23423729 27.88542857 27.752249049 28.1025 28.2415 28.37 28.05366667 28.6565 28.7755 28.898 29.013];
v0=[0.8;0.001;0.1;]; % Initial values
% Creating a matrix
zsol = [];
v1sol = [];
v2sol = [];
v3sol = [];
for i=1:numel(rho) - 1
rho0=17.1;
g=9.8;
zspan = [z(i) z(i+1)];
Nsquare = (g/rho0)*(rho(i+1)-rho(i))/(z(i+1)-z(i));
options=odeset('RelTol',1e-4,'AbsTol',1e-6,'Events',@events);
[t,v] = ode45(@(t,v)rhs(t,v,Nsquare), zspan, v0,options);
v0 = [v(end,1) ; v(end,2) ; v(end,3)];
zsol = [zsol;t];
v1sol = [v1sol;v(:,1)];
v2sol = [v2sol;v(:,2)];
v3sol = [v3sol;v(:,3)];
end
plot(v1sol,zsol,'r',v2sol,zsol,'g',v3sol,zsol,'b')
xlabel('Width,b and vertical velocity,w')
ylabel('Height, z')
grid on ;
function parameters=rhs(~,v,Nsquare)
alpha=0.116;
db= 4*alpha*v(2).^2-v(1).*v(3)./2*v(2).^2;
dw= v(1).*v(3)-4*alpha*v(2).^4+v(1).*v(2).^2.*v(3)./2*v(1).*v(2).^3;
dgmark= (-2*Nsquare*v(1).*v(2)^4-v(1).*v(3)^2+4*alpha*v(2).^4.*v(3)-v(1).*v(2).^2.*v(3).^2-8*alpha*v(2).^3.*v(3)+2*v(3).^2.*v(1).*v(2))./2*v(1).*v(2).^4;
parameters=[db;dw;dgmark];
end
  2 个评论
Dereje
Dereje 2017-12-30
yes your are right. I posted the original code. here is the event codes I used in the program.
options=odeset('RelTol',1e-4,'AbsTol',1e-6,'Events',@events);
[t,v] = ode45(@(t,v)rhs(t,v,Nsquare), zspan, v0,options);
%
%
function [value, isterminal, direction] = events(t,v)
value = v(2);
isterminal = 1;
direction = 0;
end

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Jan
Jan 2017-12-30
编辑：Jan 2017-12-30
But your works sufficiently already. Insert this test:
[t,v] = ode45(@(t,v)rhs(t,v,Nsquare), zspan, v0,options);
if t(end) ~= zspan
disp('Event function has stopped the integration');
end
You will see, that the event function works as expected. At i=14 the 0 is reached by v(2). So what does "couldn't get result" mean?
Note that you have some other troubles:
Warning: Failure at t=1.802503e+02. Unable to meet integration tolerances
without reducing the step size below the smallest value allowed
(4.547474e-13) at time t.
It seems like there is a pole. Then the integration cannot work reliably.
  5 个评论
显示 3更早的评论
Dereje
Dereje 2018-1-3
Also When I say It can not stop what I meant was: let's say v(2) is zero at z=80. Though in the result graph there is values for z greater than 80 (In principle I would have to see results only for z 0 to 80).

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by