# ode45 for sinc function

1 view (last 30 days)
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!

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
Ok, i will try to implement another control inputs. Thanks a lot!

### Categories

Find more on Programming in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by