Hello all. I am trying to solve this equation for theta(t) and plot the result:
It describes the rotation of a generator shaft (angular position of the shaft is theta, angular velocity is dtheta/dt, etc) constrained by a spring, subject to dynamic friction mju, and driven by a sinusoidal forcing function. Note EGR is a single constant, not 3 separate constants. The solution (theta(t) aka position) and its derivative (dtheta(t)/dt aka angular velocity) should be sinusoidal as well (unless some of my parameters are off).
I am trying to use ode45, which should be fine as this function is continuous, however it seems to run forever. I tried a similar problem that I ran across in this post and it completed in seconds, so I started with that code.
I've noticed by playing around that ode45 will solve it in seconds if I get rid of the 0.5*m*r^2 in the denominator (after rearranging the above for the 2nd derivative of theta, see code below), but then of course my result is off. I'm not sure why this makes such a difference in the amount of time the solver takes, because this amounts to just dividing by a constant...anyone know?
I've tried other solvers like ode15s, ode23s, 23t, 23tb, with mixed and not very reasonable results depending on how I set the parameters and tolerances in options.
I've checked and rechecked my equation/code to make sure my input is correct.
Anyway, I'm stuck at this point, so if anyone has any advice on how to solve this and plot the result, I would really appreciate it. Here's my code if you want to try it out:
m = 0.001;
r = 0.005;
EGR = 1;
Q = (0.003531/513.3)/2;
u = 0.1;
k = 30.7;
taustall = 0.003531;
xi = 0.02;
B = 6/(2*pi);
C = (3*pi)/2;
D = k.*xi;
A = (2.2-D);
kinematic = @(t,y) ([ y(2,1); ((((((A/2)*sin((B*t)+C))+(A/2)+D)*r))-((sign(y(2,1))*(EGR*((-Q)*abs(y(2,1))+taustall))+(k*(r^2)*y(1,1))+(k*xi*r)+(sign(y(2,1))*u*y(2,1)))))/((0.5*m*(r^2)))]);
[t y] = ode45(kinematic, [0:0.1:20], [0;0]);
plot(t, y(:,1), '-b', 'LineWidth',2)
plot(t, y(:,2), '-r', 'LineWidth',1)
Update: I just tried some other solvers, and ode23t with 'RelTol' set to 1e-12 gives an interesting, but not-quite-believable plot, up until about t=5, where it goes completely insane. Could be a clue, though.
Update: Perhaps the sign() is causing weird behavior with ode45 due to some discontinuity or problem with step sizes as I've seen Jan and others mention. I thought this would not be a problem because my function is continuous (right?). I'm not sure how to implement the Events function or the necessary loop to restart solving at points of discontinuity with new initial conditions, and it's kind of irrelevant until I can get some manner of solution in the first place. Can anyone shed light on this or give me an example (beyond the help documentation) of how I might implement such a loop in this case?