Plot ODE45 data in a surface

2 次查看（过去 30 天）
Mikel
Mikel 2022-7-27
回答： Star Strider 2022-7-27
Ran in:
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 个评论
Torsten
Torsten 2022-7-27
编辑：Torsten 2022-7-27
Explain in more detail what exactly you are trying to plot over [-4,4] x [-4,4].

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Star Strider
Star Strider 2022-7-27
Ran in:
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.
.

类别

在 Help Center 和 File Exchange 中查找有关 Surface and Mesh Plots 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by