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 views (last 30 days)
Show older comments
%%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 Comments
Walter Roberson
on 30 Dec 2017
You did not post code for events and the MATLAB function named events would not be relevant to this matter ?
Accepted Answer
Jan
on 30 Dec 2017
Edited: Jan
on 30 Dec 2017
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 Comments
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!