special criteria for ode solver
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
The function alpha is defined from -1 to 1 because of the arcussinus
How can I make this work so that the ode solver changes the formula if 2.52 is reached?
Any ideas?
clear
Parameter
syms phi(t) m c g alpha(t) CB_x CB_z BA_x BA_z v_x v_z x0
m= 100;
copt= 3558.7;
g= 9.81;
x0=1.342007;
k=0.75;
phi_d=diff(phi);
phi_dd= diff(phi,t,2);
if phi(t) < 2.52
alpha(t)=asin((1.7-2.1*cos(phi(t)))/3.6);
elseif 2.52<=phi(t)<=2.53
alpha(t)=0
elseif phi(t)>=2.53
asin((1.7+2.1*cos(phi(t)))/3.6);
end
CB_x= -2.1*sin(phi);
CB_z= 2.1*cos(phi);
BA_x= 1.8*cos(alpha(t));
BA_z= 1.8*sin(alpha(t));
v_x=diff(CB_x+BA_x);
v_z=diff(CB_z+BA_z);
alpha_d=diff(alpha(t),t);
v=sqrt(v_x^2+v_z^2);
% v=sqrt(2.1^2*phi_d^2+1.8^2*alpha_d^2-2*2.1*1.8*phi_d*alpha_d*sin(phi-alpha))
r1=((CB_x+BA_x)^2+(CB_z+BA_z)^2)^0.5;
x=sqrt(2.1^2+0.6^2-2*2.1*0.6*cos(phi(t)))-0.3105829;
U=1/2*copt*k*(x-x0)^2;
W=-m*g*1.8*sin(alpha(t));
V=U+W;
T_trans=0.5*m*v^2;
T_rot=1/3*m*3.6^2*alpha_d^2;
T=T_trans+T_rot;
L_1=diff(diff(T,phi_d),t);
L_2=diff(T,phi(t));
L_3=diff(V,phi(t));
% L_1=diff(diff(T,phi_d),phi)*phi_d+diff(diff(T,phi_d),phi_d)*phi_dd;
phi_dd_test=(diff(diff(T,phi_d),phi)*phi_d-L_3+L_2)/diff(diff(T,phi_d),phi_d);
F=L_1-L_2+L_3;
[Y,S] = odeToVectorField(F==0)
tspan=[0 30];
y0=[deg2rad(35.9) 0];
M = matlabFunction(Y,'vars', {'t','Y'});
[t,Y] = ode45(@(t, Y) M(t,Y),tspan,y0);
plot(t, Y, 'linewidth', 1.5)
xlim([0.00 10.00])
ylim([-4.00 4.00])
zlim([-1.00 1.00])
legend(["Winkel","Winkelgeschwindigkeit"])
xlabel("Time")
3 commentaires
Torsten
le 19 Juin 2022
@Patrick Nowohradsky comment moved here:
Yes u are right there should be an alpha=
Sam Chak
le 20 Juin 2022
Can you confirm if your alpha looks like this? Only the real parts are plotted from 0 to .
Note that the principal branch of (arcsine) is . The branch cuts into the complex plane at this interval, or approximately .
Do you really want to make at this interval ?
Réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!