Using Runge Kutta to solve second-order differential equations with two degrees of freedom
10 views (last 30 days)
Show older comments
Hi everyone, I am beginner in Matlab Programming,I have been trying to solve the differential equation system given below, but I am unable to establish a workflow. Any assistance or related work would be greatly appreciated. Thanks in advance.
In the following system of differential equations, x and xb are structural responses, omga is frequency ratio, tau is time, and other parameters are constants. The requirement is that we need to provide the relationship between the structural response of the system and the frequency ratio omga.
matlab program:
function [T,X,dX] = ODE_RK4( Hfun,t,h,x0 )
if nargin < 4
error('The initial value must be given');
end
if isstr(Hfun)
eval(['Hfun = @',Hfun,';']);
end
n = length(t);
if n == 1
T = 0:h:t;
elseif n == 2
T = t(1):h:t(2);
else
T = t;
end
T = T(:);
N = length(T);
x0 = x0(:);
x0 = x0';
m = length(x0);
X = zeros(N,m);
dX = zeros(N,m);
X(1,:) = x0;
for k = 2:N
h = T(k) - T(k-1);
K1 = Hfun( T(k-1) , X(k-1,:)' );
K2 = Hfun( T(k-1)+h/2 , X(k-1,:)'+h*K1/2 );
K3 = Hfun( T(k-1)+h/2 , X(k-1,:)'+h*K2/2 );
K4 = Hfun( T(k-1)+h , X(k-1,:)'+h*K3 );
X(k,:) = X(k-1,:)' + (h/6) * ( K1 + 2*K2 + 2*K3 + K4 );
dX(k-1,:) = (1/6) * ( K1 + 2*K2 + 2*K3 + K4 );
end
dX(N,:) = Hfun( T(N),X(N,:) );
if nargout == 0
plot(T,X)
end
clc;
clear;
lambda_b = 0.01;
lambda_ns = -0.01;
lambda = 0.5;
zeta_b = 0.025;
chi = 0.0001;
alpha = 0.2;
zeta = 0.05;
z = 0.06;
omega_values = linspace(0.05, 1, 1000);
tau = linspace(0, 2*pi, 1000);
dt = tau(2) - tau(1);
X0 = [0; 0; 0; 0]; % Initial value
x_max = zeros(1, length(omega_values));
for i = 1:length(omega_values)
omega = omega_values(i);
f = @(t, X) [X(2);
-(2*zeta*omega*X(2) + 4*(2*alpha^2*X(1)^3*lambda+1/sqrt(alpha)) + lambda_b*(X(1) - X(3)) + z*cos(t)) / omega^2;
X(4);
-1/(chi*omega^2)*(2*zeta_b*omega*X(4) + lambda_ns*X(3) - lambda_b*(X(1) - X(3)))];
[T,X]= ODE_RK4(f, tau, dt,X0);
x_max(i) = max(abs(X(i, :)));
end
figure;
plot(omega_values, x_max);
xlabel('omega');
ylabel('x_{max}');
0 Comments
Answers (2)
Torsten
on 30 Sep 2024
Edited: Torsten
on 30 Sep 2024
I don't understand from your code if you want to perform a parameter study for 1000 different values of omega or if you want to define omega as a function of t.
It looks like a parametric sweep, but this line needs some explanations:
x_max(i) = max(abs(X(i, :)));
tspan = [0 2*pi];
y0 = [0; 0; 0; 0]; % Initial value
Omega = @(t)0.05+t*(1-0.05)/(2*pi);
[T,Y] = ode45(@(t,y)fun(t,y,Omega),tspan,y0);
plot(T,Y)
function dydt = fun(t,y,Omega)
% y(1) = x_b, y(2) = x_b_dot, y(3) = x, y(4) = x_dot
omega = Omega(t);
lambda_b = 0.01;
lambda_ns = -0.01;
lambda = 0.5;
zeta_b = 0.025;
chi = 0.0001;
alpha = 0.2;
zeta = 0.05;
z = 0.06;
dydt = zeros(4,1);
dydt(1) = y(2);
dydt(2) = (-2*zeta_b*omega*y(2)-lambda_ns*y(1) + lambda_b*(y(3)-y(1)))/(chi*omega^2);
dydt(3) = y(4);
dydt(4) = (-2*zeta*omega*y(4)-4*(2*alpha^2*y(3)^3*lambda+1/sqrt(alpha))-chi*omega^2*dydt(2)-...
2*zeta_b*omega*y(2)-lambda_ns*y(1)-z*cos(t))/omega^2;
end
0 Comments
Alan Stevens
on 30 Sep 2024
Like this?
omega = 0.05:0.05:20;
tspan = [0 2*pi];
xmax = zeros(1,numel(omega));
for i = 1:numel(omega)
X0 = [0,0,0,0];
[t, X] = ode45(@(t,X) fn(t,X,omega(i)),tspan,X0);
x = X(:,1);
xmax(i) = max(abs(x));
end
plot(omega,xmax),grid
xlabel('omega'), ylabel('xmax')
function dXdt = fn(t, X, omega)
zeta = 0.05;
alpha = 0.2;
lambda = 0.5;
lambda_b = 0.01;
chi = 0.0001;
zeta_b = 0.025;
lambda_ns = -0.01;
z = 0.06;
x = X(1); v = X(2); xb = X(3); vb = X(4);
dXdt = [v;
(-z*cos(t) - lambda_b*(x - xb)...
- 4*(2*alpha^2*x^3*lambda+1/sqrt(alpha))...
-2*zeta*omega*v)/omega^2;
vb;
(lambda_b*(x-xb) - lambda_ns*xb...
- 2*zeta_b*omega*vb)/(chi*omega^2)];
end
0 Comments
See Also
Categories
Find more on Spectral Estimation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!