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

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}');

Answer 1

[Torsten]

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

Answer 2

[Alan Stevens]

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