Using Runge Kutta to solve second-order differential equations with two degrees of freedom

10 views (last 30 days)
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}');

Answers (2)

Torsten
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

Alan Stevens
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

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!