No solution for my script/seemingly infinite calculation time

I need to calculate and graph solutions for a grain on a sifting plate, and I have this script, and it works for most of my values, but for some reason when i use these values it just keeps on calculating, and matlab starts using more of my RAM

close all

clear all

clc

R=0.025;                % [m]

Omega=4;                % [rad/s]

Mu=0.5;                 % [dimensionless]

g=9.81;                 % [m/s^2]

tspan=[0 100];               % time interval

ystart=[0 0.0000001 0 0];    % beginning conditions

eq=@(t,y) [y(2) ; - (Mu*g*y(2))/(sqrt( (y(2)^2) + (y(4)^2) )) + R*(Omega^2)*(cos(Omega*t)) ; y(4) ;  - (Mu*g*y(4))/(sqrt( (y(2)^2) + (y(4)^2) )) + R*(Omega^2)*(sin(Omega*t))  ];

[t,y]=ode45(eq,tspan,ystart);

figure(1)
plot(y(:,1), y(:,3), 'g');     % this plots the position

xlabel('positie along the x-axis');

ylabel('positie along the y-axis');

title('movement of grain combination 2');

axis equal;

[Walter Roberson]

Generally speaking that would tend to suggest that it is approaching a singularity. For example y2 and y4 could both be going to 0

Answer 1

[Walter Roberson]

Your y(2) is oscillating badly compared to its absolute values, so the ode is having to take a lot of small steps to do the integration.

Have a look at

- (Mu*g*y(2))/(sqrt( (y(2)^2) + (y(4)^2) )) + R*(Omega^2)*(cos(Omega*t))

R*(Omega^2)*(cos(Omega*t)) is close to being constant near t = 0, a value roughly 0.4

In - (Mu*g*y(2))/(sqrt( (y(2)^2) + (y(4)^2) )), the sign of y(4) is not relevant because y(4) only appears ^2. Whether y(2) is positive or negative, (sqrt( (y(2)^2) + (y(4)^2) )) is going to be positive provided that at least one of y(2) or y(4) is non-zero. Mu and g are positive. With y(2) being positive, (Mu*g*y(2))/(sqrt( (y(2)^2) + (y(4)^2) )) would be positive. But then there is the leading "-" which would take that sub-expression negative.

So we are now adding a negative expression to R*(Omega^2)*(cos(Omega*t)) . What happens if the result is negative? Well, if it is, then next time through y(2) would be negative, and (Mu*g*y(2))/(sqrt( (y(2)^2) + (y(4)^2) )) would be negative, and with the leading "-" that would get you positive, and added to the positive R*(Omega^2)*(cos(Omega*t)) the result would be positive. And next time around that might lead you to negative... and around and around you go.

Because this switch between negative and positive would be happening on every step, ode45 has to take small steps.