ode45 for sinc function

5 views (last 30 days)
Muhammad
Muhammad on 20 Mar 2023
Commented: Muhammad on 20 Mar 2023
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!

Accepted Answer

Sam Chak
Sam Chak on 20 Mar 2023
Edited: Sam Chak on 20 Mar 2023
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
  3 Comments
Sam Chak
Sam Chak on 20 Mar 2023
Edited: Sam Chak on 20 Mar 2023
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
Muhammad on 20 Mar 2023
Ok, i will try to implement another control inputs. Thanks a lot!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!