In the following program replace the function of duhamel3 with a Newmark beta method function, the results of the program must be the same as using the duhamel3 function. Please provide the code of the Newmark beta method in matlab and a brief explanation.

%---------------------- Program VIBRA1.m ----------------------------% % % % Program to calculate the response of a viscously underdamped single% % DOF system to an arbitrary force F(t). using using recursive % % equation. % % % %------------------ Last updated: November 12 -2017 -----------------% clear all; close all %----------- System properties --------------------------------------% m = 0.85; % mass of the system (N.m.s^2/rad) c = 0.9; % damping coefficient (N.m.s/rad) k = 140.0; % stiffness coefficient (N.m/rad) u0 = 0.0; % initial angular displacement (rad) v0 = 0.0; % initial angular velocity (rad/s) dt = 0.0002; % time interval of the external moment (sec) F0 = 4; % amplitude of the moment (N.m) tf = 10; % final time(sec) %-------- Read and plot force time history ------------------------% t = (0:dt:tf)'; nt = length(t); for ii =1:nt ft(ii)=0; end for ii = 1:10000 ft(ii) = F0*sin(3.14/2*t(ii)); end figure; plot( t, ft ); grid on; xlabel('Time [sec]'); ylabel('Moment [N.m]') %---------- Calculate the natural frequency, natural period and damping ratio -----------% wn = sqrt( k / m ); Tj = 2*pi ./ wn; zj = c / (2* m * wn); wd = wn *sqrt(1-zj^2); disp(['*** The natural frequency [rad/s] is :',num2str(wn)]); disp(['*** The natural period [sec] is :',num2str(Tj)]); disp(['*** The damping ratio is :',num2str(zj)]); disp(['*** The damped natural frequency [rad/s] is :',num2str(wd)]); %---------- Calculate the modal and physical response ----------- [u , v , a] = duhamel3(wn,zj,m,dt,u0,v0,ft,t); % Duhamel variation linear force %---------- Plot the response time histories -------------------- disp(['*** The maximum rotation [rad.] is :',num2str(max(u))]); figure; plot( t,u, 'MarkerSize',10 ); grid on legend('Angular Displacement') title('Angular Displacement time history '); ylabel('Angular displacement [rad]'); xlabel('Time [sec]') figure; plot( t,v, 'MarkerSize',10 ); grid on legend('Angular Velocity') title('Angular velocity time history '); ylabel('Angular velocity [rad/s]'); xlabel('Time [sec]') figure; plot( t,a, 'MarkerSize',10 ); grid on legend('Angular Acceleration') title('Angular acceleration time history '); ylabel('Angular acceleration [rad/s^2]'); xlabel('Time [sec]')

------------------------------------------------------------------------------------------------- function [u,v,a] = duhamel3( wn, z, m, dt, u0, v0, f, t ) % Program to calculate the response of a linear viscously damped % SDOF system subjected to an excitation f(t) with arbitrary time % variation sampled at constant intervals of duration 'dt'. % % The excitation is assumed to be constant in each time interval. % % d^2(q(t))/dt^2 + 2*z*w*d(q(t))/dt + (w^2)*q(t) = f(t) % % d^2(q(t))/dt^2 = f(t) - 2*z*w*d(q(t))/dt - (w^2)*q(t) % % Input data: % ---------- % w = Natural frequency [rad/sec]. % z = Damping ratio for underdamped case. % m = Mass of oscillator [consistent units]. % dt = Time interval [sec]. % u0 = Initial displacement [Length]. % v0 = Initial velocity [Length/sec]. % f = Vector with the sampled excitation [consistent units]. % t = Vector with discrete time [sec]. % % Output data: % ----------- % u = Displacement time history. % v = Velocity time history. % a = Acceleration time history. % Initialize: nt = length(t); Z = zeros(3,nt); % Initial conditions: Z(1,1) = u0; Z(2,1) = v0; Z(3,1) = ( f(1) / m ) - (wn^2)*u0 - 2*z*wn*v0; % Calculate coefficients: wd = wn * sqrt(1-z^2); ex = exp(-z*wn*dt); co = cos(wd*dt); si = sin(wd*dt); a11 = ( co + z/sqrt(1-z^2) * si ) * ex; a12 = 1 / wd * si * ex; a21 = -wn^2 * a12; a22 = ( co - z/sqrt(1-z^2) * si ) * ex; SS = sqrt(1-z^2); kas = 1/(m*dt*SS*wn^3); b11 = kas*(2*z*SS+((1-2*z^2-z*wn*dt)*si-(2*z*SS+wd*dt)*co)*ex); b12 = kas*(wd*dt-2*z*SS+(-(1-2*z^2)*si+2*z*SS*co)*ex); b21 = kas*(-wd+(wn*(z+wn*dt)*si+wd*co)*ex); b22 = kas*(wd-(wn*z*si+wd*co)*ex); A1 = [a11 , a12 ; a21 , a22] ; B1 = [b11 , b12 ; b21 , b22] ; % Compute response: for i = 2 : nt-1 ff = [f(i-1) ; f(i)]; Z(1:2,i) = A1 * Z(1:2,i-1) + B1 * ff; Z(3,i) = ( f(i) / m ) - (wn^2)*Z(1,i) - 2*z*wn*Z(2,i); end u = Z(1,:); v = Z(2,:); a = Z(3,:); return