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