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;
1 个评论
Walter Roberson
2017-3-26
Generally speaking that would tend to suggest that it is approaching a singularity. For example y2 and y4 could both be going to 0
回答（1 个）
Walter Roberson
2017-3-26
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.
类别
在 帮助中心 和 File Exchange 中查找有关 Numeric Solvers 的更多信息
产品
2017-3-26
2017-3-26
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
