Plotting different datasets from looped ode45 using subplots.
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Frederik Bjerregaard
am 9 Dez. 2021
Kommentiert: Star Strider
am 9 Dez. 2021
Hi. I have a script where I'm solving and ode for different values of a parameter. I want to plot each of these cases on a seperate subplot as plotting them all on the same plot is messy. So far i have managed to get the result for the first iteration onto every subplot using:
clc; clear;close all;
% Defining parameters and initial conditions
f = 0.35+[0.1:0.1:0.4];
r = 0.3;
zeta = 0.25;
eps = -1/6;
Theta0 = [0.30, 0];
tspan = [0 100];
% Solver options
opts = odeset('RelTol',1e-6,'AbsTol',[1e-9 1e-9]);
% Solver
for i = 1:length(f)
[t,theta] = ode45(@(t,theta) odefcn(t,theta,f(i),r,zeta), tspan,Theta0);
eqn=@(A) f(i)^2-(-r^2*A+A+(3*eps*A^3)/4)^2-(2*zeta*r*A)^2;
Amp(:,i)=fzero(eqn,0.1);
thetaHB=Amp(i).*sin(r.*t);
theta_HB=Amp(i).*r.*cos(r.*t);
subplot(2,2,1);
plot(theta(:,1),theta(:,2),thetaHB,theta_HB);
pbaspect([1 1 1])
title('f=0.35+0.1');
xlabel ('\theta');
ylabel('d\theta/dt');
legend('Numerical','Harmonic Balance');
subplot(2,2,2);
plot(theta(:,1),theta(:,2),thetaHB,theta_HB);
pbaspect([2 2 1])
title('f=0.35+0.2');
xlabel ('\theta');
ylabel('d\theta/dt');
legend('Numerical','Harmonic Balance');
subplot(2,2,3);
plot(theta(:,1),theta(:,2),thetaHB,theta_HB);
pbaspect([1 1 1])
title('f=0.35+0.3');
xlabel ('\theta');
ylabel('d\theta/dt');
legend('Numerical','Harmonic Balance');
subplot(2,2,4);
plot(theta(:,1),theta(:,2),thetaHB,theta_HB);
pbaspect([1 1 1])
title('f=0.35+0.4');
xlabel ('\theta');
ylabel('d\theta/dt');
legend('Numerical','Harmonic Balance');
sgtitle('Phase Planes r=0.3','FontSize',20)
end
The underlying function for odefcn is
function dthetadt = odefcn(t,theta,f,r,zeta)
dthetadt = zeros(2,1);
dthetadt(1)=theta(2);
dthetadt(2)=f.*sin(r.*t)-2.*zeta.*theta(2)-sin(theta(1));
end
0 Kommentare
Akzeptierte Antwort
Star Strider
am 9 Dez. 2021
I assume that the intended result is to have different plots representing the function with different values of ‘f’ in each subplot. The code needs to be tweaked dlightly to get there —
% Defining parameters and initial conditions
f0 = 0.35;
fv = [0.1:0.1:0.4];
r = 0.3;
zeta = 0.25;
eps = -1/6;
Theta0 = [0.30, 0];
tspan = [0 100];
% Solver options
opts = odeset('RelTol',1e-6,'AbsTol',[1e-9 1e-9]);
% Solver
for i = 1:numel(fv)
f = f0 + fv(i); % Create New 'f' For Each Iteration
[t,theta] = ode45(@(t,theta) odefcn(t,theta,f,r,zeta), tspan,Theta0);
eqn=@(A) f^2-(-r^2*A+A+(3*eps*A^3)/4)^2-(2*zeta*r*A)^2;
Amp(:,i)=fzero(eqn,0.1);
thetaHB=Amp(i).*sin(r.*t);
theta_HB=Amp(i).*r.*cos(r.*t);
subplot(3,2,i); % Iterate Through 'subplot' Array
plot(theta(:,1),theta(:,2),thetaHB,theta_HB);
pbaspect([1 1 1])
title(sprintf('f=0.35+%.1f',fv(i)));
xlabel ('\theta');
ylabel('d\theta/dt');
hl = legend('Numerical','Harmonic Balance');
hl.Visible = 'off';
% subplot(2,2,2);
% plot(theta(:,1),theta(:,2),thetaHB,theta_HB);
% pbaspect([2 2 1])
% title('f=0.35+0.2');
% xlabel ('\theta');
% ylabel('d\theta/dt');
% legend('Numerical','Harmonic Balance');
%
% subplot(2,2,3);
% plot(theta(:,1),theta(:,2),thetaHB,theta_HB);
% pbaspect([1 1 1])
% title('f=0.35+0.3');
% xlabel ('\theta');
% ylabel('d\theta/dt');
% legend('Numerical','Harmonic Balance');
%
% subplot(2,2,4);
% plot(theta(:,1),theta(:,2),thetaHB,theta_HB);
% pbaspect([1 1 1])
% title('f=0.35+0.4');
% xlabel ('\theta');
% ylabel('d\theta/dt');
% legend('Numerical','Harmonic Balance');
% sgtitle('Phase Planes r=0.3','FontSize',20)
end
subplot(3,2,[5 6])
ax = gca;
pos = ax.Position;
hl.Position = pos + [0.1 0.15 -0.2 -0.2];
hl.Visible = 'on';
ax.Visible = 'off';
sgtitle('Phase Planes r=0.3','FontSize',20)
function dthetadt = odefcn(t,theta,f,r,zeta)
dthetadt = zeros(2,1);
dthetadt(1)=theta(2);
dthetadt(2)=f.*sin(r.*t)-2.*zeta.*theta(2)-sin(theta(1));
end
I re-positioned the legend so that it would not cover the plots, since it was common to all of them.
.
2 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Subplots finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!