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

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

[DGM]

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.

Answer 1

[Jan]

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.

[Ugen Kezang]

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

Answer 2

[Davide Masiello]

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

[Jan]

Fine! Thanks.