フィルターのクリア

Plot ODE45 data in a surface

1 回表示 (過去 30 日間)
Mikel
Mikel 2022 年 7 月 27 日
回答済み: Star Strider 2022 年 7 月 27 日
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 日
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 ExchangeScatter Plots についてさらに検索

製品


リリース

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by