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.
2 次查看(过去 30 天)
显示 更早的评论
%%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 个评论
Walter Roberson
2017-12-30
You did not post code for events and the MATLAB function named events would not be relevant to this matter ?
采纳的回答
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!