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

2 次查看（过去 30 天）
Ugen Kezang
Ugen Kezang 2022-4-29
编辑： Davide Masiello 2022-4-29
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 个评论
DGM
DGM 2022-4-29
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Jan
Jan 2022-4-29
编辑：Jan 2022-4-29
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 个评论
Ugen Kezang
Ugen Kezang 2022-4-29
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

请先登录，再进行评论。

Davide Masiello
Davide Masiello 2022-4-29
编辑：Davide Masiello 2022-4-29
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by