Can someone help me understand why I cant get a graph with an alpha of 0.001

1 view (last 30 days)
Ugen Kezang
Ugen Kezang on 29 Apr 2022
Edited: Davide Masiello on 29 Apr 2022
I'm unable to get a plot for an alpha of 0.001 as I'm trying to plot behaviour of theta with alpha at 0.1, 0.01 and 0.001
clear;
clc;
close all;
Question_number = 4.2;
if Question_number == 4.2
[t,x] = ode15s(@dc_motor_42,[0 0.5],[1;1;0]);
plot(t,x(:,1),'-o');
end
function dxdt = dc_motor_42(t,x)
J = 3.2284E-6;
b = 3.5077E-6;
K = 0.0274;
R = 4;
L = 2.75E-6;
Kc = 0; % Set Kc to 5.
u_dist = 0; % Set different voltage input
alpha = 0.1;
if x(3) > alpha
dxdt = [0*x(1) 1*x(2) 0*x(3); ...
0*x(1) -b/J*x(2) K/J*(x(3)-alpha); ...
-Kc/L*x(1) -K/L*x(2) -R/L*x(3)]+[0;0;u_dist/L];
elseif x(3) < -alpha
dxdt = [0*x(1) 1*x(2) 0*x(3); ...
0*x(1) -b/J*x(2) K/J*(x(3)+alpha); ...
-Kc/L*x(1) -K/L*x(2) -R/L*x(3)]+[0;0;u_dist/L];
else
dxdt = [0 1 0; ...
0 -b/J 0; ...
-Kc/L -K/L -R/L]*x + [0;0;u_dist/L];
end
end
  1 Comment
DGM
DGM on 29 Apr 2022
Not sure, but I see that your function returns dxdt as 3x1 in the third case due to the multiplication with x, but the other two cases return a 3x3 matrix.
It looks like it only starts entering those cases when alpha = 0.001 or so.

Sign in to comment.

Answers (2)

Jan
Jan on 29 Apr 2022
Edited: Jan on 29 Apr 2022
With alpha = 0.1 only the 3rd case of the if-branches occurs. Then dxdt is replied as [3 x 1] vector.
With alpha = 0.001 the first and second branch can occur also. Then dxdt is a [3 x 3] matrix. In the 3rd case you multiply a [3x3] matrix by a [3x1] vector, but in the first two cases it is an elementwise multiplication between a [3x3] matrix and a [1x3] vector:
if x(3) > alpha
dxdt = [0*x(1) 1*x(2) 0*x(3); ...
0*x(1) -b/J*x(2) K/J*(x(3)-alpha); ...
-Kc/L*x(1) -K/L*x(2) -R/L*x(3)] + [0;0;u_dist/L];
elseif x(3) < -alpha
dxdt = [0*x(1) 1*x(2) 0*x(3); ...
0*x(1) -b/J*x(2) K/J*(x(3)+alpha); ...
-Kc/L*x(1) -K/L*x(2) -R/L*x(3)] + [0;0;u_dist/L];
else
dxdt = [0 1 0; ...
0 -b/J 0; ...
-Kc/L -K/L -R/L] * x + [0;0;u_dist/L];
end
You want this instead:
if x(3) > alpha
dxdt = [0*x(1) + 1*x(2) + 0*x(3); ...
0*x(1) + -b/J*x(2) + K/J*(x(3)-alpha); ...
-Kc/L*x(1) + -K/L*x(2) + -R/L*x(3)] + [0;0;u_dist/L];
elseif x(3) < -alpha
dxdt = [0*x(1) + 1*x(2) + 0*x(3); ...
0*x(1) + -b/J*x(2) + K/J*(x(3)+alpha); ...
-Kc/L*x(1) + -K/L*x(2) + -R/L*x(3)] + [0;0;u_dist/L];
else
dxdt = [0 1 0; ...
0 -b/J 0; ...
-Kc/L -K/L -R/L] * x + [0;0;u_dist/L];
end
But a simplification would be easier:
if x(3) > alpha
v = K / J * (x(3) - alpha);
elseif x(3) < -alpha
v = K / J * (x(3) + alpha);
else
v = 0;
end
dxdt = [0 1 0; ...
0 -b/J v; ...
-Kc/L -K/L -R/L] * x + [0; 0; u_dist/L];
Finally consider, that your function to be integrated is not smooth. Matlab's ODE integrators are designed to handle smooth functions only, so this is not a reliable method anymore. The numerically correct version is to stop the integration, when the limits are reached and to restart it. Your method might produce a result, but it can be dominated by rounding errors.
This foum contains hundreds of questions concerning the integration of non-smooth functions. Unfortunately Matlab's documentation contains an example also, which uses interp1 to create a time dependent parameter. This is junk from the viewpoint of numerical maths. Do not use such methods for a scientific publication.
  1 Comment
Ugen Kezang
Ugen Kezang on 29 Apr 2022
Thank you very much for this. I realised the matrices were not of compatible sizes later on, but couldn't figure out how to manipulate the matrix into the required size

Sign in to comment.


Davide Masiello
Davide Masiello on 29 Apr 2022
Edited: Davide Masiello on 29 Apr 2022
Following @Jan comprehensive answer, I think a code like the one suggested below might be a good solution to your problem
clear, clc
alpha = 0.001;
J = 3.2284E-6;
b = 3.5077E-6;
K = 0.0274;
R = 4;
L = 2.75E-6;
Kc = 0;
u_dist = 0;
x1 = 1;
x2 = 1;
x3 = 0;
t0 = 0;
tf = 0.5;
region = 2;
idx = 0;
while t0 < tf
idx = idx+1;
if region == 2 && x3 == alpha
A = [ 0 1 0 ;...
0 -b/J K/J ;...
-Kc/L -K/L -R/L ];
c = [ 0; -K/J*alpha; u_dist/L];
elseif region == 2 && x3 == -alpha
A = [ 0 1 0 ;...
0 -b/J K/J ;...
-Kc/L -K/L -R/L ];
c = [ 0; K/J*alpha; u_dist/L];
else
A = [ 0 1 0 ;...
0 -b/J 0 ;...
-Kc/L -K/L -R/L ];
c = [ 0; 0; u_dist/L];
end
options = odeset('Event',@(t,x)eventLocator(t,x,alpha));
sol(idx) = ode15s(@(t,x)dc_motor_42(t,x,A,c),[t0 tf],[x1;x2;x3],options);
t0 = sol(idx).x(end);
x1 = sol(idx).y(1,end);
x2 = sol(idx).y(2,end);
x3 = sol(idx).y(3,end);
if x3 > alpha
region = 3;
elseif x3 < -alpha
region = 1;
else
region = 2;
end
end
plot([sol(:).x],[sol(:).y],'-o');
legend('x1','x2','x3')
function dxdt = dc_motor_42(t,x,A,c)
dxdt = A*x+c;
end
function [position,isterminal,direction] = eventLocator(t,x,alpha)
position = abs(x(3))-alpha; % The value that we want to be zero
isterminal = 1; % Halt integration
direction = 0; % The zero can be approached from either direction
end

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!