special criteria for ode solver
1 view (last 30 days)
Show older comments
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 Comments
Torsten
on 19 Jun 2022
@Patrick Nowohradsky comment moved here:
Yes u are right there should be an alpha=
Sam Chak
on 20 Jun 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 ?
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!