Plot ODE45 data in a surface

3 views (last 30 days)
Mikel
Mikel on 27 Jul 2022
Answered: Star Strider on 27 Jul 2022
Hello,
I have a 2 order system solved with ode 45 which I want to somehow plot it into the plane that I created using meshgrid. I did some trials using contour, plot3 and surf but I didnt had success. Can you help me?
this is my code:
%% INIT
clear variables;
close all;
clc;
%% variable initialization
% time
time=150;
t = (0:0.1:time);
%coefs
a=-0.3;
b=-1.1;
%initial conditions
x1init= 3;
x2init= 0;
% vertex
d1=-4;
d2=4;
%% SOLVE SYSTEM AND CREATE THE PLANE
[t,x] = ode45(@(t,x) systemfunc(a,b,x), t, [x1init x2init]);
x1=x(:,1)';
x2=x(:,2)';
% working area
[X,Y]=meshgrid(d1:2:d2);
% dimensions of the system
X1=X(1,:);
X2=Y(:,1);
%% PLOT
%plot system trayectory
plot(x1,x2,'k', 'LineWidth',2)
set(gca, 'FontSize', 12)
hold on
plot(x1(1,1),x2(1,1),'bo','LineWidth',2) % starting point
plot(x1(1,end),x2(1,end),'ks','LineWidth',2) % ending point
% plot position and speed of the system
figure
subplot(2,1,1)
plot(t,x1,'k', 'LineWidth',2)
set(gca, 'FontSize', 12)
hold on
xlabel('time(s)','fontweight','bold','FontSize',12)
ylabel('position','fontweight','bold','FontSize',12)
subplot(2,1,2)
plot(t,x2,'k', 'LineWidth',2)
set(gca, 'FontSize', 12)
xlabel('time(s)','fontweight','bold','FontSize',12)
ylabel('speed','fontweight','bold','FontSize',12)
title('function evolution')
%plot surface generated by the function
Z=@(X,Y)(a*X+b*Y);
figure
s=surf(X,Y,Z(X,Y));
set(gca, 'FontSize', 12)
s.FaceColor = [0 0.7 1];
%view(0,90)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FUNCTIONS
function dxdt=systemfunc(a,b,x)
dxdt=zeros(2,1);
dxdt(1) = x(2);
dxdt(2) = a*x(1)+b*x(2);
end
  1 Comment
Torsten
Torsten on 27 Jul 2022
Edited: Torsten on 27 Jul 2022
Explain in more detail what exactly you are trying to plot over [-4,4] x [-4,4].

Sign in to comment.

Answers (1)

Star Strider
Star Strider on 27 Jul 2022
One option for creating a surface from two vectors is to multiply the transpose of one with the other so that the first vector is a column vector and the second vector is a row vector. If both are positive, then an accurate representation would involve taking the element-wise square root of the resulting matrix, however if there are any negative values, that reuslts in a complex matrix (as is the problem here), so that would need to be dealt with by either not calculating the square root, or by taking the absolute values or plotting the real and imaginary parts separately —
%% INIT
clear variables;
close all;
clc;
%% variable initialization
% time
time=150;
t = (0:0.1:time);
%coefs
a=-0.3;
b=-1.1;
%initial conditions
x1init= 3;
x2init= 0;
% vertex
d1=-4;
d2=4;
%% SOLVE SYSTEM AND CREATE THE PLANE
[t,x] = ode45(@(t,x) systemfunc(a,b,x), t, [x1init x2init]);
x1=x(:,1)';
x2=x(:,2)';
% working area
[X,Y]=meshgrid(d1:2:d2);
% dimensions of the system
X1=X(1,:);
X2=Y(:,1);
%% PLOT
%plot system trayectory
plot(x1,x2,'k', 'LineWidth',2)
set(gca, 'FontSize', 12)
hold on
plot(x1(1,1),x2(1,1),'bo','LineWidth',2) % starting point
plot(x1(1,end),x2(1,end),'ks','LineWidth',2) % ending point
% plot position and speed of the system
figure
subplot(2,1,1)
plot(t,x1,'k', 'LineWidth',2)
set(gca, 'FontSize', 12)
hold on
xlabel('time(s)','fontweight','bold','FontSize',12)
ylabel('position','fontweight','bold','FontSize',12)
subplot(2,1,2)
plot(t,x2,'k', 'LineWidth',2)
set(gca, 'FontSize', 12)
xlabel('time(s)','fontweight','bold','FontSize',12)
ylabel('speed','fontweight','bold','FontSize',12)
title('function evolution')
%plot surface generated by the function
z = (x1(:)*x2);
% % Z=@(X,Y)(a*X+b*Y);
figure
s = surfc(t, t, abs(z), 'EdgeColor','none');
xlim([0 15])
ylim([0 15])
% s=surf(X,Y,Z(X,Y));
% % set(gca, 'FontSize', 12)
% % s.FaceColor = [0 0.7 1];
% view(0,90)
% %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FUNCTIONS
function dxdt=systemfunc(a,b,x)
dxdt=zeros(2,1);
dxdt(1) = x(2);
dxdt(2) = a*x(1)+b*x(2);
end
I commented-out the view call. Restore it to see this as a sort of contour plot.
I am not certain what result you want.
.

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!