6 views (last 30 days)
Show older comments
Muhammad on 20 Mar 2023
-
-
Link
Direct link to this question
https://support.mathworks.com/matlabcentral/answers/1931740-ode45-for-sinc-function
Commented: Muhammad on 20 Mar 2023
Accepted Answer: Sam Chak
Open in MATLAB Online
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!
0 Comments Show -2 older commentsHide -2 older comments
Show -2 older commentsHide -2 older comments
Sign in to comment.
Sign in to answer this question.
Accepted Answer
Sam Chak on 20 Mar 2023
Edited: Sam Chak on 20 Mar 2023
Open in MATLAB Online
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
3 Comments Show 1 older commentHide 1 older comment
Show 1 older commentHide 1 older comment
Muhammad on 20 Mar 2023
Direct link to this comment
https://support.mathworks.com/matlabcentral/answers/1931740-ode45-for-sinc-function#comment_2669235
Edited: Muhammad on 20 Mar 2023
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 on 20 Mar 2023
Direct link to this comment
https://support.mathworks.com/matlabcentral/answers/1931740-ode45-for-sinc-function#comment_2669515
Edited: Sam Chak on 20 Mar 2023
Open in MATLAB Online
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 on 20 Mar 2023
Direct link to this comment
https://support.mathworks.com/matlabcentral/answers/1931740-ode45-for-sinc-function#comment_2669695
Ok, i will try to implement another control inputs. Thanks a lot!
Sign in to comment.
More Answers (0)
Sign in to answer this question.
See Also
Categories
Control SystemsControl System ToolboxLinear AnalysisTime and Frequency Domain AnalysisTime-Domain Analysis
Find more on Time-Domain Analysis in Help Center and File Exchange
Tags
- ode45
- sinc
- control systems
- variable structure control
- nonlinear control systems
Products
- Control System Toolbox
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)
Contact your local office
