Solution of two order differential equations; Fourth order Runge-Kutta method iteration

Hi,i solute this two order differential equations in fourth order Runge-Kutta method iteration. However, my calculation results are different from the amplitude in the literature.

This is the two order differential equations:

is the acceleration, velocity and displacement response of the structure. All three quantities are unknown, and the former is the derivative of the latter.

is the relative acceleration, velocity and displacement between the damper and the structure. All three quantities are unknown, and the former is the derivative of the latter.

Ms'、meq are quality, Ks、keq are stiffness, Cs、ceq are damping. But ceq is not a fixed value, but a quantity that varies with xr.

F(t) is the external load and is a sin function with amplitude of F0=11.7 and frequency of fz. The frequency is unknown.

The result in the paper is the result of normalization. The x-coordinate is beta = fz / 3.51; The y-coordinate is Xs.max / (F0 / Ks).

In this case ,there are three main variables, namely the two displacements Xs and xr, and the frequency of the external load fz.

Firstly, I simplified this second order differential equation in matrix form:

Secondly, i used ‘Maple' to simplify these two equations:

So, i get 4 state equations.

Thirdly, i used the fourth order Runge-Kutta method iteration to get the result..

But my results don't match the literature, and I don't know why.

I want to know if there is any error in my iterative calculation operation. If not, maybe I set wrong parameters when I read the literature

Mss = 4062.9;  %Mss is the Ms'                                       
Ks = 49656;
Cs = 14;
meq = 77.6;
keq = 913;
F0 = 11.7;                                         


t_end=100;  %The total time in the time domain
n=10000;    %Total number of samples in the time domain
h=t_end/n;  %Step in the time domain
time=linspace(0,100,n); %Generating time sequence
m=100;      %Total number of samples in the frequency domain
beta=linspace(0.7,1.3,m); %Generate frequency domain sequence

for i=1:m
    betai=beta(i);   %Single extraction of data in the frequency domain
    fz=3.51*betai;  %External excitation frequency
    y0=zeros(4,1);  %The initial value of the equation of state [Xs,Vs,Xt,Vt]
    y=y0;           %Assign the initial value to the equation of state once per beta cycle
    collection(:,1)=y;  %Collect initial state
    for j=1:n           %Fourth order Runge Kutta method
        t=time(j);
        K1=g3(t,y,fz);
        K2=g3(t+h/2,y+h/2*K1,fz);
        K3=g3(t+h/2,y+h/2*K2,fz);
        K4=g3(t+h,y+h*K3,fz);
        y=y+h/6*(K1+2*K2+2*K3+K4);
        collection(:,j)=y;  %Record the response
    end
    xs=collection(1,:);
    xr=collection(3,:);
    maxd(i)=max(abs(xs))/(F0/Ks);
end
figure()
plot(beta,maxd)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy=g3(t,y,fz)
ceq = 1006.4 * y(3);
meq = 77.6;
keq = 913;
Mss = 4062.9; 
Ks = 49656;
Cs = 14;
F = 11.7 * sin(fz * t);    %Excitation form:sin(t)

%Import equation of state
dy = zeros(4,1);
dy(1) = y(2);
dy(2) = (ceq*y(4)+keq*y(3)-Cs*y(2)-Ks*y(1)+F)/Mss;
dy(3) = y(4);
dy(4) = -(y(4)*ceq*meq+y(4)*ceq*Mss+y(3)*keq*meq+y(3)*keq*Mss-Cs*y(2)*meq-Ks*y(1)*meq+F*meq)/(Mss*meq);
%
end

Answer 1

[Torsten]

Running your problem with ODE45 gives your Runge-Kutta results. So the problem seems to be that your parameters are different from those in the publication.

Mss = 4062.9;  %Mss is the Ms'                                       
Ks = 49656;
Cs = 14;
meq = 77.6;
keq = 913;
F0 = 11.7;                                         
t_end=100;  %The total time in the time domain
n=10000;    %Total number of samples in the time domain
h=t_end/n;  %Step in the time domain
time=linspace(0,100,n); %Generating time sequence
m=100;      %Total number of samples in the frequency domain
beta=linspace(0.7,1.3,m); %Generate frequency domain sequence
for i=1:m
    betai=beta(i);   %Single extraction of data in the frequency domain
    fz=3.51*betai;  %External excitation frequency
    y0=zeros(4,1);  %The initial value of the equation of state [Xs,Vs,Xt,Vt]
    y=y0;           %Assign the initial value to the equation of state once per beta cycle
    collection(:,1)=y;  %Collect initial state
    [T,SOL] = ode45(@(t,y)g3(t,y,fz),time,y0);
    xs = SOL(:,1);
    maxd(i)=max(abs(xs))/(F0/Ks);
end
figure()
plot(beta,maxd)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy=g3(t,y,fz)
ceq = 1006.4 * y(3);
meq = 77.6;
keq = 913;
Mss = 4062.9; 
Ks = 49656;
Cs = 14;
F = 11.7 * sin(fz * t);    %Excitation form:sin(t)
%Import equation of state
dy = zeros(4,1);
dy(1) = y(2);
dy(2) = (ceq*y(4)+keq*y(3)-Cs*y(2)-Ks*y(1)+F)/Mss;
dy(3) = y(4);
%dy(4) = -(y(4)*ceq*meq+y(4)*ceq*Mss+y(3)*keq*meq+y(3)*keq*Mss-Cs*y(2)*meq-Ks*y(1)*meq+F*meq)/(Mss*meq);
dy(4) = (-meq*dy(2)-ceq*y(4)-keq*y(3))/meq;
%
end

[Yijing LU]

Thanks a lot! Your answer is very useful to me.