PLEASE HELP!!!!!!! How do i modify this code so that i use the Runge-Kutta metho

Question

PLEASE HELP!!!!!!! How do i modify this code so that i use the Runge-Kutta method so that the code can be used to solve the solution of the force vibration of a single DOF system with viscous damping......................

% clear command window, variables, and opened figure windows

clc, clear, close all

% set up the plot window

ScreenSize = get(0,'ScreenSize');
FigVibration = figure('Name','Vibration', ...
                      'NumberTitle','off', ...
                      'OuterPosition',[0.1*ScreenSize(3), 0.1*ScreenSize(4) 0.75*ScreenSize(3) 0.75*ScreenSize(4)]);
% input parameters
              
Mass = 100;
Stiffness = 4e5;
Zeta = 0.1;
InitialDisplacement = 0;
InitialVelocity = 0;
FMax = 100;
Frequency = sqrt(Stiffness/Mass);

% parameters

NaturalFrequency = sqrt(Stiffness/Mass);
Period = 2*pi/Frequency;
FM = FMax/Mass;
CC = 2*Zeta*NaturalFrequency;
CK = NaturalFrequency^2;
DTime = Period/(2^14);
XMax = (FMax/Stiffness)/(2*Zeta);

% initial conditions

V(1) = InitialVelocity;
X(1) = InitialDisplacement;
Time(1) = 0;
i = 1;
j = 1;

plot(Time/Period,X/XMax,'-or','linewidth',2);
title('The Forced Vibration')
xlabel('time [sec]')
ylabel('X')
hold on

% update subsequent results using the Euler method

while Time(i) < 20*Period
    Time(i+1) = Time(i) + DTime;
    V(i+1) = V(i) + (FM*sin(Frequency*Time(i))-CC*V(i)-CK*X(i))*DTime;
    X(i+1) = X(i) + V(i)*DTime;
    i = i+1;
    if Time(i)-Time(j)>0.05*Period
        plot(Time(j:i)/Period,X(j:i)/XMax,'-or','linewidth',2);
        j = i;
    end
    drawnow
end

% plot the approximate harmonic solution against the accurate result


% just some lines for me to see conveniently if the results are accurate

line([0 20], [1 1])
line([0 20], [0 0])
line([0 20], [-1 -1])

Explanation / Answer

RUNGE-KUTTA

Forward Euler

function [u,v] = feuler(a,b,h,ya)
u= (a:h:b)';
n=length(u);
v=ya(:)';
for i=1:n-1,
v(i+1,:) = v(i,:) + h*Ode_function(u(i),v(i,:));
end

Heun

function [u,v] = heun(a,b,h,ya)
u= (a:h:b)';
n=length(u);
v=ya(:)';
for i=1:n-1,
w= v(i,:) + h*Ode_function(u(i),v(i,:));
v(i+1,:) = v(i,:)+ (.5)*h*(Ode_function(u(i),v(i,:))+Ode_function(u(i+1),w));
end

4th Order RK

function [u,v] = rk4(a,b,h,ya)
u= (a:h:b)';
n=length(u);
v=ya(:)';
for i=1:n-1,
u0=u(i); h2=h/2;
v0=v(i,:); f0=Ode_function(u0,v0);
v1=v0+h2*f0; f1=Ode_function(u0+h2,v1);
v2=v0+h2*f1; f2=Ode_function(u0+h2,v2);
v3=v0+h*f2; f3=Ode_function(u0+h,v3);
v(i+1,:) = v0+(h/6)*(f0+2*f1+2*f2+f3);
end

4th Order RKF

Related Questions

Navigate

• Browse (All)
• Browse P
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.