HW4 solutions - (1 Convert the equation given to a system...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: (1) Convert the equation given to a system of differential equations that can be solved using the methods we have learned and using MATLAB. y= 1 2 1 2 1 2 2 2 1 (2) Write a function that contains this system of differential equations. function deriv = springDerHw4_solutions(t,y) % First we unpack the variables in the y vector into something without % parentheses y1 = y(1); y2 = y(2); % Establish our constant values c = 5; % N*s/m k = 20; % N/m m = 20; % kg deriv = zeros(2,1); % Now we simply implement the derivative equations. deriv(1) = y2; deriv(2) = (‐c*y2 ‐ k*y1)/m; (3) Euler’s method clear all % Define our knowns: t0 = 0; tf = 15; x0 = 1; v0 = 0; h = .1; % We can define our t vector because we have a set step size and defined % limit of integration tVec = [t0:h:tf]; % Allocate our other vectors x = zeros(size(tVec)); v = zeros(size(tVec)); %Initialize x(1) = x0; v(1) = v0; % Perform Euler integration % NOTE: this loop goes to length(tVec) MINUS ONE because the loop defines % step (i+1); so when i gets to length ‐ 1, the next defined value is the % end of the array. If we went the whole length of the tVec, then we'd % have a longer concentration and temperature vector than time vector. for i=1:length(tVec) ‐ 1 derivs = springDerHw4_solutions(tVec(i),[x(i) v(i)]); dxdt = derivs(1); dvdt = derivs(2); x(i+1) = x(i) + h*dxdt; v(i+1) = v(i) + h*dvdt; end (4) Compute MATLAB’s solutions with ode45 % Now we do ODE45 [rkTime, rkY] = ode45(@springDerHw4_solutions,[t0 tf],[x0 v0]); rkX = rkY(:,1); rkV = rkY(:,2); (5) Plot solutions of Euler’s method and ode45 in the same figure, and save the figure as a PNG file. % Now we plot it all up figure(1) subplot(2,1,1) plot(tVec,x,'b‐') hold on plot(rkTime,rkX,'k‐') xlabel('time (s)') ylabel('position (m)') title('Integrated mass position') legend('Euler','Runge‐Kutta') hold off subplot(2,1,2) plot(tVec,v,'b‐') hold on plot(rkTime,rkV,'k‐') xlabel('time(s)') ylabel('velocity (m/s)') legend('Euler','Runge‐Kutta') title('Integrated mass velocity') hold off suptitle('Comparing Euler and R‐K for a damped spring‐mass system') The above picture is saved as PNG file. (6) See if h=0.1 is small enough to be accurate. (Answer as comment lines) % Part 5: h = 0.01 is definitely the best and reproduces what ode45 can do. % h = 0.1 is OK... it shows damping, but has some trouble getting the % amplitude of the oscillations right. % As we go to h = 0.25, we see that we have bigger problems, as the % oscillations don't seem to damp, but stay constant instead, so the % magnitude of the error gets worse as time goes on. % At h = 0.5, things are clearly getting much worse, as now the system is % diverging from what it should be, and the amplitude and velocity are % increasing with oscillations. % When integrated to a final time of 150s, we see that the solution % explodes and is clearly very far away from what it should be (the % amplitude of the position and velocity are on the order of 10^7. ...
View Full Document

This note was uploaded on 05/01/2011 for the course CHBE 2120 taught by Professor Gallivan during the Spring '07 term at Georgia Tech.

Ask a homework question - tutors are online