ode45 for sinc function

Hey guys, i want to ask if ode45 can be used to solve state equation that has a sinc function in it.

My state equation is:

function dstatesdt = midterm2(t,state)

    x1 = state(1);
    x2 = state(2);

    z1 = 1/(5^2 - x1^2);
    z2 = 1/(5^2 - x2^2);

%     u = - (x1*(z1/z2 + 1) + z1*x1^2*sinc(x2)/z2 + x2^2 + x2);x
    u = - (x1*(z1/z2 + 1) + 2*z1*x1^2*sinc(x2)/z2 + x2^2 + x2);
%     u = - (x1*(z1/z2 + 1) + z1*x1^2*sinc(x2)/z2 + x2^2 + x2 + x2/z1);

    x1dot = x2 + x1*sin(x2);
    x2dot = x2^2 + x1 + u;

    dstatesdt = [x1dot;x2dot];

And my function to solve and plot the ODE is:

function sim_ode(dstatesdt,tspan,bx1,dx1,bx2,dx2,arg)
    
    for x1 = bx1(1):dx1:bx1(2)
        for x2 = bx2(1):dx2:bx2(2)
            xi = [x1,x2] %show initial point
            [to,stateo] = ode45(dstatesdt,tspan,xi);
            
            x1o = stateo(:,1);
            x2o = stateo(:,2);

            % if any(x1o >= 5) | any(x1o <= -5) | any(x2o >= 5) | any(x2o <= -5)
            %     fprintf('Overbound for Initial State [%s,%d] \n',x1,x2)
            % end

            if arg == 1
                plot(to,x1o)
            elseif arg == 2
                plot(to,x2o)
            else
                plot(x1o,x2o)
%                 plot(x1,x2,'bo')
%                 plot(x1o(end),x2o(end),'ks')
            end
            drawnow
        end
    end
end

And i run this code on main script:

f = @midterm2;

figure(1)
hold on
title('Time Series of x_1')
xlabel('t (s)')
ylabel('x_1 (t)')
grid()
sim_ode(f,[0 10],[-4.5 4.5],0.5,[-4.5 4.5],0.5,1)
hold off

at first inital point [-4.5 -4.5] the codes didn't plot anything and didn't show any error and just kept running (around 5 min) so i stop it and this showed up in command prompt:

xi =

   -4.5000   -4.5000

Operation terminated by user during ode45


In sim_ode (line 6)
            [to,stateo] = ode45(dstatesdt,tspan,xi);

In main_sim (line 15)
sim_ode(f,[0 10],[-4.5 4.5],0.5,[-4.5 4.5],0.5,1)

Any help would be appriciated!

Thanks!

Answer 1

[Sam Chak]

Hi @Muhammad

Although your system is highly nonlinear, the ode45 can evaluate the sinc() function. However, your choice of initial values causes the trajectory to diverge to , where the singularity (division-by-zero) occurs in this term .

[x, y] = meshgrid(-4.9:0.35:4.9);
z1 = 1./(5^2 - x.^2);
z2 = 1./(5^2 - y.^2);
w = - (x.*(z1./z2 + 1) + 2*z1.*(x.^2).*sinc(y)./z2 + y.^2 + y);
u = y + x.*sin(y);
v = y.^2 + x + w;
l = streamslice(x, y, u, v);
axis tight

Choosing another set of initial values, for example, causes the trajectory to converge to the origin.

tspan  = linspace(0, 10, 101);
x0     = [-4.5 4.5];
[t, x] = ode45(@midterm2, tspan, x0);

plot(t, x, 'linewidth', 1.5), grid on
xlabel('t'), ylabel('\bf{x}'), 
legend('x_{1}', 'x_{2}')

function dstatesdt = midterm2(t, state)

    x1 = state(1);
    x2 = state(2);

    z1 = 1/(5^2 - x1^2);
    z2 = 1/(5^2 - x2^2);

%     u = - (x1*(z1/z2 + 1) + z1*x1^2*sinc(x2)/z2 + x2^2 + x2);x
    u = - (x1*(z1/z2 + 1) + 2*z1*x1^2*sinc(x2)/z2 + x2^2 + x2);
%     u = - (x1*(z1/z2 + 1) + z1*x1^2*sinc(x2)/z2 + x2^2 + x2 + x2/z1);

    x1dot = x2 + x1*sin(x2);
    x2dot = x2^2 + x1 + u;

    dstatesdt = [x1dot; x2dot];
end

[Muhammad]

is this means that i choose the wrong control input? for the context, i have to find a control input to make the system asymptotically stable at the origin (0,0) which constrained by for any initial condition in that boundary

[Sam Chak]

Hi @Muhammad

Your originally query on the sinc() function is cleared up. However, I'm unfamiliar with your design of u. Your system is highly nonlinear and due to this trigonometric term , there can be multiple equilibrium points. So, I suggest that a Variable Structure Control approach to be implemented. You can test different initial values.

If you find the additional solution and suggested control laws are helpful, please consider accepting ✔ and voting 👍 the Answer. Thanks a bunch! 🙏

[x, y] = meshgrid(-4.9:0.35:4.9);
% w = (- 2*(x.*sin(y) + y) - 2*x - y.^2 - (x.*sin(y) + y).*sin(y) - x.*cos(y).*(y.^2 + x))./(x.*cos(y) + 1);
w = - 5*tanh((x + y)/0.1) - (y.^2 + x);
u = y + x.*sin(y);
v = y.^2 + x + w;
l = streamslice(x, y, u, v); hold on
plot(0, 0, 'ro', 'linewidth', 2, 'MarkerSize', 10), hold on
plot(-pi/2, pi/2, 'ro', 'linewidth', 2, 'MarkerSize', 10), hold on
plot(3*pi/2, -3*pi/2, 'ro', 'linewidth', 2, 'MarkerSize', 10), hold off
axis tight

tspan  = linspace(0, 30, 3001);
% x0     = [-5 -5];     % Test 1
% x0     = [-5  5];     % Test 2
x0     = [ 5 -5];       % Test 3
% x0     = [ 5  5];     % Test 4
[t, x] = ode45(@midterm2, tspan, x0);

plot(t, x, 'linewidth', 1.5), grid on
xlabel('t'), ylabel('\bf{x}'), 
legend('x_{1}', 'x_{2}')

function dstatesdt = midterm2(t, state)

    x = state(1);
    y = state(2);
    
    % Variable Structure Control 
    if (x >= -pi/2) && (y <= pi/2)
        u = (- 2*(x*sin(y) + y) - 2*x - y^2 - (x*sin(y) + y)*sin(y) - x*cos(y)*(y^2 + x))/(x*cos(y) + 1);
    elseif (x <= 3*pi/2) && (y >= -3*pi/2)
        u = (- 2*(x*sin(y) + y) - 2*x - y^2 - (x*sin(y) + y)*sin(y) - x*cos(y)*(y^2 + x))/(x*cos(y) + 1);
    else
        u = - 5*tanh((x + y)/0.1) - (y^2 + x);
    end

    x1dot = y + x*sin(y);
    x2dot = y^2 + x + u;

    dstatesdt = [x1dot; x2dot];
end

[Muhammad]

Ok, i will try to implement another control inputs. Thanks a lot!