ode45 giving incorrect result

I'm trying to control a pendulum to be stable when inverted and my phase portrait seems to show that it's stable there (the origin). (X axis is angle Y axis is angular velocity. But when i use ode45 on the function that produced that phase portrait it gives weird results.

Here's the script, and the function pendups that the script calls:

close all;

%Phase portrait bounds and density
lx=-pi;
dx=0.1;
ux=pi;

ly=-5;
dy=0.1;
uy=5;

g=9.81;
l=1;
r=pi/2;


%Phase portrait creation
[X1,X2]=meshgrid(lx:dx:ux,ly:dy:uy);
F1=X1;
F2=X2;
for i=1:length(X1(:,1))
    for j=1:length(X1(1,:))
        x=[X1(i,j);X2(i,j);l;r];
        f=pendups(0,x);
        F1(i,j)=f(1);
        F2(i,j)=f(2);
    end
end

%Plot phase portrait
figure;
quiver(X1,X2,F1,F2);

%Plot specific solution (should be tangent to phase portait)

x=[-1,2,l,r];
[xt,x_traj]=ode45('pendups', [0 10], x')
hold on;
plot(x_traj(:,1),x_traj(:,2))
xlim([lx,ux]);
ylim([ly,uy]);
hold off;

Now the function:

function f=pendups(t,x)
x=x';
g=9.81;
f=[0,0,0,0];
f(1)=x(2);
%control state (close to inverted with low enough velocity):
if(sqrt(x(1)^2+x(2)^2)<x(4))
    f(2)=g/x(3)*sin(x(1))-(2*g/x(3)*x(1)+g/x(3)/2*x(2))*cos(x(1));
%uncontrolled state:
else
    f(2)=g/x(3)*sin(x(1));
end
f(3)=x(3);
f(4)=x(4);
f=f';
end

Here's the result, the phase portrait looks correct but the specific solution obviously is not tangent to it. This happens with other initial conditions as well, some better some worse.

Is this an issue with the function being discontinuous? If so is there a way to fix it?

Answer 1

[Sam Chak]

Hi @Joseph Billips

There were errors in your code when you attempt to describe the dynamics of the pendulum, which is a 2nd-order system, using 4 states. I guess that l = 1 and r = pi/2 are constants. Check if the response should look like this:

% Phase portrait bounds and density
lx = -pi;
ux =  pi;
dx =  (ux - lx)/20;
ly = -5;
uy =  5;
dy =  (uy - ly)/20;
g  =  9.81;
l  =  1;
r  =  pi/2;

% Phase portrait creation
[X1, X2] = meshgrid(lx:dx:ux, ly:dy:uy);
F1 = X1;
F2 = X2;
for i = 1:length(X1(:, 1))
    for j = 1:length(X1(1, :))
        x = [X1(i, j); X2(i, j)];
        f = pendups(0, x);
        F1(i,j) = f(1);
        F2(i,j) = f(2);
    end
end

% Plot phase portrait
figure;
quiver(X1, X2, F1, F2);

% Plot specific solution (should be tangent to phase portait)
% x0 = [-1 2 l r];
x0 = [-1 2];
[xt, x_traj] = ode45(@pendups, [0 5], x0);
hold on;
% plot(xt, x_traj(:,1))
plot(x_traj(:,1), x_traj(:,2), 'linewidth', 1.5)
xlim([lx, ux]);
ylim([ly, uy]);
hold off;
xlabel('x_{1}(t)')
ylabel('x_{2}(t)')

function f = pendups(t, x)
    x  =  x';
    lx = -pi;
    dx =  0.1;
    ux =  pi;
    ly = -5;
    dy =  0.1;
    uy =  5;
    g  =  9.81;
    l  =  1;
    r  =  pi/2;
    % f  =  [0, 0, 0, 0];
    f  =  [0, 0];
    
    f(1) = x(2);
    
    % control state (close to inverted with low enough velocity):
    if(sqrt(x(1)^2 + x(2)^2) < r)
        f(2) = g/l*sin(x(1)) - (2*g/l*x(1) + ((g/l)/2)*x(2))*cos(x(1));
    % uncontrolled state:
    else
        f(2) = g/l*sin(x(1));
    end
    
%     % control state (close to inverted with low enough velocity):
%     if(sqrt(x(1)^2 + x(2)^2) < x(4))
%         f(2) = g/x(3)*sin(x(1)) - (2*g/x(3)*x(1) + ((g/x(3))/2)*x(2))*cos(x(1));
%     % uncontrolled state:
%     else
%         f(2) = g/x(3)*sin(x(1));
%     end
    
%     f(3) = x(3);
%     
%     f(4) = x(4);
    
    f = f';
end