# ode45 for sinc function

조회 수: 1(최근 30일)
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 2023년 3월 20일
편집: Sam Chak 2023년 3월 20일
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표시 이전 댓글 수: 2숨기기 이전 댓글 수: 2
Ok, i will try to implement another control inputs. Thanks a lot!

댓글을 달려면 로그인하십시오.

### 범주

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!