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

6 views (last 30 days)
Yijing LU on 6 Jul 2022
Commented: Yijing LU on 7 Jul 2022
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

Torsten on 6 Jul 2022
Edited: Torsten on 6 Jul 2022
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 on 7 Jul 2022