How to save dydt(3) equation computed inside ODE?
3 ビュー (過去 30 日間)
古いコメントを表示
Hi, i have to use an ode45 with a function with some if and elseif conditions inside it. The idea is that the ODE system [dydt(1) ; dydt(2) ; dydt(3)] does change depending on the solutions y(1), y(2) and y(3) that the ode45 gives each time istant of integration. The code works the way i writed it below, but now i have to plot the dydt(3) array and i dont know how. I thought to use plot( t(1:end-1), diff(y(:,3)) ), and it also work but it gives a very bad solution (too different from the real one). I also have seen the solutions proposed here:
but it doesnt work for me and i dont know why. There is a way to save the dydt(3) that result from my function? Thanks in advance.
MATLAB Version: 9.11.0.1837725 (R2021b) Update 2 .
[t,y] = ode45(@(t,y) odefcn(t,y), [0 T_f], zeros(1,3)) ;
function dydt = odefcn(t,y)
global omega_0 N_Orbits T T_f J_0 m L sigma zeta K_d K k_theta tau_theta ...
theta_REF theta_REF_dot T_M SAT
dydt = zeros(3,1) ;
if y(3) == 0
dydt(1) = y(2);
dydt(2) = cos(t) - (k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot )) ;
dydt(3) = k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot ) ;
elseif abs(y(3)) > SAT
if sign( y(3) ) == sign( k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot ) )
dydt(1) = y(2);
dydt(2) = exp(-t) ;
dydt(3) = 0 ;
elseif sign( y(3) ) ~= sign( k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot ) )
dydt(1) = y(2);
dydt(2) = cos(t) - (k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot )) ;
dydt(3) = k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot ) ;
end
elseif abs(y(3)) <= SAT
dydt(1) = y(2);
dydt(2) = cos(t) - (k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot )) ;
dydt(3) = k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot ) ;
end
end
4 件のコメント
Walter Roberson
2022 年 5 月 27 日
Yes there is, but it is almost always the wrong thing to do. ode45 evaluates the function at a variety of locations in order to predict how the function evolves, and then it tests the prediction against a different location. If the test and predictions are too different then the step is rejected and a smaller step is used; upon success a larger step will be used. It follows that except for systems that evolve as a polynomial of low degree, that ode45 will increase the step size until there is a prediction failure, so failed steps are normal. Therefore the saved history is pretty much guaranteed to include rejected locations. This is not what you want to plot.
Jan
2022 年 5 月 27 日
@Alessandro Castello: Yes, this was a typo in my former answer. It is fixed now:
[t, y] = ode45(@fcn, [0, 2*pi], [0, 0]);
[~, p] = fcn(t.', y.');
function [dy, p] = fcn(t, y)
p = sin(y(2, :)); % <- The wanted parameter
dy = [2 .* t; p ./ y(1, :)];
end
This method does not have the problem Walter mentions. It is evaluated for the acepted steps only. Using the OutputFcn is valid also.
回答 (1 件)
Torsten
2022 年 5 月 27 日
[t,y] = ode45(@(t,y) odefcn(t,y), [0 T_f], zeros(1,3)) ;
dydt = zeros(size(y));
for i = 1:numel(t)
dydt(i,:) = odefcn(t(i),y(i,:));
end
plot(t,dydt(:,3))
0 件のコメント
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!