
How can I plot a graph with a time dependent ode45?
    4 visualizzazioni (ultimi 30 giorni)
  
       Mostra commenti meno recenti
    
    Donghun Lee
 il 16 Mar 2020
  
    
    
    
    
    Commentato: Ameer Hamza
      
      
 il 16 Mar 2020
            clc
clear all
A = 0.05; %Excitation Amplitude
l_r = 2; %Wave length of the road
v = 45; %Speed(m/s)
P = l_r/v; %Period
Om = 1/P*2*pi; %Forcing Frequency
%Om = %0.07m,2m,45m/s
k_l = 26400; %Linear stiffness
m = 483; %Mass
d = -0.1; %Stretching condition
l = 0.5; %Length of the spring
k_s = -(k_l*(l-d))/(4*d); %Spring stiffness
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%% 
f = @(t,x) [ x(2); ...
    -(4*k_s*(x(1)-A*sin(Om*t))*(sqrt((l-d)^2 + (x(1)-A*sin(Om*t))^2) - l))/ ...
    (m*(sqrt((l-d)^2 + (x(1)-A*sin(Om*t))^2))) ];
T = 50;
x0 = [0,0];
[t,x] = ode45(f,[0,T],x0);
Response_amp = max(x(:,1)) - min(x(:,1));
plot(t,x(:,1))
xlabel('Time (s)')
ylabel('Amplitude (m)')
title('When d=-0.1', 'fontsize', 20)
set(gca,'FontSize',15)
Hi, all. This is my ode45 code. If I chagne A (Excitation Amplitude), Response_amp will have a different value. I wish to plot this relationship as A vs Response_amp.
0 Commenti
Risposta accettata
  Ameer Hamza
      
      
 il 16 Mar 2020
        
      Modificato: Ameer Hamza
      
      
 il 16 Mar 2020
  
      You can visualize such relation with 3D plots
l_r = 2; %Wave length of the road
v = 45; %Speed(m/s)
P = l_r/v; %Period
Om = 1/P*2*pi; %Forcing Frequency
%Om = %0.07m,2m,45m/s
k_l = 26400; %Linear stiffness
m = 483; %Mass
d = -0.1; %Stretching condition
l = 0.5; %Length of the spring
k_s = -(k_l*(l-d))/(4*d); %Spring stiffness
f_n = sqrt(k_l/m)/(2*pi); %Natural frequency
%%
fig = figure();
ax = axes();
hold(ax);
view([-53 33]);
grid on
A_array = 0.05:0.05:0.3; %Excitation Amplitude
T = 15;
x0 = [0,0];
for i=1:numel(A_array)
    A = A_array(i);
    f = @(t,x) [ x(2); ...
        -(4*k_s*(x(1)-A*sin(Om*t))*(sqrt((l-d)^2 + (x(1)-A*sin(Om*t))^2) - l))/ ...
        (m*(sqrt((l-d)^2 + (x(1)-A*sin(Om*t))^2))) ];
    [t, x] = ode45(f,[0,T],x0);
    Response_amp = max(x(:,1)) - min(x(:,1));
    plot3(t, A*ones(size(t)), x(:,1), 'LineWidth', 1);
end
xlabel('Time (s)')
ylabel('A (m)')
zlabel('Amplitude (m)')
title('When d=-0.1')

2 Commenti
Più risposte (0)
Vedere anche
Categorie
				Scopri di più su Numerical Integration and Differential Equations in Help Center e File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!