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.
- Spring '07