special criteria for ode solver

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")



 

[John D'Errico]

What do you think this line does?

asin((1.7+2.1*cos(phi(t)))/3.6);

Hint: NOTHING. It produces a result, then dumps it into the bit bucket.

I would note that in the other branches of that loop, you actually assign something to alpha.

Should i get into questions like, can you assign a SPECIFIC value to a symbolic function at a point? NO. So lines like this are meaningless anyway.

alpha(t)=asin((1.7-2.1*cos(phi(t)))/3.6);

You might want to learn about the piecewise utility in the symbolic toolbox.

[Torsten]

@Patrick Nowohradsky comment moved here:

Yes u are right there should be an alpha=

[Sam Chak]

Hi @Patrick Nowohradsky

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 ?