%------------------------------------------------------------------------ %%% ENGINEERING MATHEMATICS III - MATH 2Z03 %%% LAB #4 %%% HIGHER ORDER ODEs: SPRING-MASS SYSTEM %------------------------------------------------------------------------ %------------------------------------------------------------------------ % Covers: % - "Advanced Engineering Mathematics" by D.G. Zill and M.R. Cullen, % Section 3.8.3 (example 6) %------------------------------------------------------------------------ %------------------------------------------------------------------------ % Instructions: % - Execute the commands given below in MATLAB. You can do this in % one of the following three ways: % (i) retype the commands one by one at the MATLAB command line (>>), % (ii) copy and paste the commands one by one to the MATLAB % command line (>>), % (iii) run all the commands at once by calling the present file % as a script, i.e. >> m2z03lab % - Analyze the output and compare it with the comments embedded in % the file below. %------------------------------------------------------------------------ clc; close all; clear all; %------------------------------------------------------------------------ % Section 3.8.3: Spring-Mass Systems: Driven Motion disp(' -----------------------------------------------------------') disp(' Section 3.8.3: Spring-Mass Systems: Driven Motion ') %%% Example: % Solve the initial-value problem % 0.2 x'' + 1.2 x' + 2x = 5cos(4t), % Initial values: % 1) x(0) = 0.5, x'(0) = 0, % 2) x(0) = 0.2, x'(0) = 1, % 3) x(0) = -1.5, x'(0) = -4, % 4) x(0) = 2, x'(0) = 2. % Exact solution (for the 1st initial conditions): % x(t) = exp(-3t)*((38/51)*cos(t)-(86/51)*sin(t))-(25/102)*cos(4t)+(50/51)*sin(4t) % After the substitution y' = u % the equation is equivalent to the system: % x' = u, % u' = -6u - 10x + 25cos(4t) % Heun_method_for_system % Initial values x01 = 0.5; u01 = 0; % 1st initial value x02 = 0.2; u02 = 1; % 2nd x03 = -1.5; u03 = -4; % 3rd x04 = 2; u04 = 2; % 4th % Initial data h = 0.05 % time step T = 2 % maximum value of t n = T/h % number of steps % Values of t, x and u at first step k = 1 (for different initial values) t(1) = 0; x1(1) = x01; u1(1) = u01; % 1st initial value x2(1) = x02; u2(1) = u02; % 2nd x3(1) = x03; u3(1) = u03; % 3rd x4(1) = x04; u4(1) = u04; % 4th % Loop to generate steps from k = 1 to k = n for k = 1:n t(k+1) = t(k) + h; % increment for t % 1st initial value x_star = x1(k) + h*u1(k); u_star = u1(k) + h*(-6*u1(k)-10*x1(k)+25*cos(4*t(k))); x1(k+1) = x1(k) + (h/2)*(u1(k)+u_star); % update for the solution x(t) u1(k+1) = u1(k) + (h/2)*(-6*u1(k)-10*x1(k)+25*cos(4*t(k))+ ... -6*u_star-10*x_star+25*cos(4*t(k+1))); % update for the solution u(t) % 2nd initial value x_star = x2(k) + h*u2(k); u_star = u2(k) + h*(-6*u2(k)-10*x2(k)+25*cos(4*t(k))); x2(k+1) = x2(k) + (h/2)*(u2(k)+u_star); % update for the solution x(t) u2(k+1) = u2(k) + (h/2)*(-6*u2(k)-10*x2(k)+25*cos(4*t(k))+ ... -6*u_star-10*x_star+25*cos(4*t(k+1))); % update for the solution u(t) % 3rd initial value x_star = x3(k) + h*u3(k); u_star = u3(k) + h*(-6*u3(k)-10*x3(k)+25*cos(4*t(k))); x3(k+1) = x3(k) + (h/2)*(u3(k)+u_star); % update for the solution x(t) u3(k+1) = u3(k) + (h/2)*(-6*u3(k)-10*x3(k)+25*cos(4*t(k))+ ... -6*u_star-10*x_star+25*cos(4*t(k+1))); % update for the solution u(t) % 4th initial value x_star = x4(k) + h*u4(k); u_star = u4(k) + h*(-6*u4(k)-10*x4(k)+25*cos(4*t(k))); x4(k+1) = x4(k) + (h/2)*(u4(k)+u_star); % update for the solution x(t) u4(k+1) = u4(k) + (h/2)*(-6*u4(k)-10*x4(k)+25*cos(4*t(k))+ ... -6*u_star-10*x_star+25*cos(4*t(k+1))); % update for the solution u(t) end % Exact solution xExact = exp(-3*t).*((38/51)*cos(t)-(86/51)*sin(t))-(25/102)*cos(4*t)+(50/51)*sin(4*t); % Plotting (exact and numerical solutions) plot(t,xExact,'-r',t,x1,'.-b'); % plotting (Figure 1) title('Forced Damped Spring-Mass System: Solution'); % figure's title, legend('exact solution','numerical approximation'); % legend xlabel('t'); ylabel('x'); % and axis label disp(' Press any key to continue ... ') pause % Plotting the phase portrait figure; % opening new window for Figure 2 plot(x1,u1,'.-b'); hold on; % 1st trajectory plot(x2,u2,'.-g'); hold on; % 2nd trajectory plot(x3,u3,'.-r'); hold on; % 3rd trajectory plot(x4,u4,'.-k'); hold off; % 4th trajectory title('Forced Damped Spring-Mass System: Phase Portrait'); % figure's title, legend(['x_{01} = ' num2str(x01) ', x\prime_{01} = ' num2str(u01)],... ['x_{02} = ' num2str(x02) ', x\prime_{02} = ' num2str(u02)],... ['x_{03} = ' num2str(x03) ', x\prime_{03} = ' num2str(u03)],... ['x_{04} = ' num2str(x04) ', x\prime_{04} = ' num2str(u04)]); % legend xlabel('x(t)'); ylabel('x\prime(t)'); % and axis label disp(' This is the end of the lab tutorial #4. ')