Download Numerical Solution of 2nd-Order ODE with Euler and Runge-Kutta - Prof. Aldo Ferri and more Assignments Computer Science in PDF only on Docsity!
ME2016, Fall 2008, Dr. Ferri Homework Assignment 12 Due Wednesday, December 3.
Problem 1. (5 points) Consider the initial-value problem
y ′ = − 2 t y^2 ; y (0) = 1 (1)
It can be shown that the exact solution to this nonlinear ode satisfying the initial condition is
y t ( ) = 1/( t^2 + 1) (2)
(a) Perform Euler integration with h = 0.1 from 0 < t < 2 sec. Make a plot of the Euler solution and the exact solution. Compute the maximum error over 0 < t < 2 sec. (b) Repeat part (a) with h = 0.01. By what factor does the maximum error decrease? Does this match your expectations? (c) Integrate the solution to (1) using the 2 nd^ -order Runge-Kutta method and h = 0.2. Plot the solution vs the exact solution and vs the Euler answer from part (a). Based on the maximum error, which is more accurate and by what factor? (d) Integrate the solution to (1) from 0 < t < 2 using Matlab’s ode45 function (Type "help ode45" for further info.) Plot the solution vs the exact solution. Compute the maximum error over 0 < t < 2 sec and compare with the errors obtained using the methodology of parts (a) and ( c).
Solution
By inspection, we see that f ( , t y ) = − 2 t y^2. For Euler integration, the first few steps are:
2 y 1 (^) = y 0 (^) + f t ( 0 (^) , y 0 ) h = 1 − 2(0)(1) (0.1) = 1 2 y 2 (^) = y 1 (^) + f t ( 1 (^) , y 1 ) h = 1 − 2(0.1)(1) (0.1) =0. 2 y 3 (^) = y 2 (^) + f t ( 2 (^) , y 2 ) h = 0.98 − 2(0.2)(0.98) (0.1) =0.
To automate the process, here is a code that generates the Euler solution, calculates the exact solution, and generates the plot.
f = @(x,y) -2xy.^2;
h = 0.1; Ya = 1; Ta = 0; tn = 0; yn = 1; while tn <=2-h;
ynp1 = yn + f(tn,yn)*h; tn = tn + h; yn = ynp1; Ta = [Ta tn]; Ya = [Ya yn];
end
y_exact_a = 1./(Ta.^2 +1);
max_error_a = max(abs(y_exact_a - Ya));
figure(1) plot(Ta,Ya,'o--',Ta,y_exact_a,'linewidth',2); xlabel('Time (sec)') ylabel('y') axis([0 2 0 1]) legend('Euler, h = 0.1','Exact') gtext(['max error = ',num2str(max_error_a)]); title('Hw12, Problem 1(a)')
You can see from the figure that the answer is fairly good, with a maximum error of 0.02632. The second figure shows a close-up of the first few steps. Figures 3 and 4 show the corresponding results for h = 0.01. We see that the computed solution is much closer to the exact solution, the maximum error reduces to 0.002416, The trend in the maximum error should follow that of the global error ; i.e., for Euler integration, the error should be of order O(h). Thus, decreasing h by a factor of 10 should reduce the error by a factor of 10. In this case, the actual reduction is by a factor of 10.9, which is very close to the prediction.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0
1
Time (sec)
y
Hw12, Problem 1(a)
max error = 0.
Euler, h = 0. Exact
Figure1. 1st^ -order Euler for h = 0.
0 0.05 0.1 0.15 0.
1
Time (sec)
y
Hw12, Problem 1(b)
Euler, h = 0. Exact
Figure 4. Close-up of Figure 3
(c) For part (c), we are to perform 2nd^ -order RK. The first few steps are as follows:
Prediction: y 1^0 = y 0 (^) + f t ( 0 (^) , y 0 ) h = 1 − 2(0)(1) (0.2)^2 =1.
Average slope:
0 2 0 0 1 1 1
f t y f t y
= ⎣^ ⎦^ = ⎣^ ⎦= −
Correction: y 1 = y 0 + φ 1 h = 1 − (0.2)(0.2) =0.
Exact: y t ( ) = 1/((0.2)^2 + 1) =0.
Prediction: y 2^0 = y 1 (^) + f t ( 1 (^) , y 1 ) h = 0.96 − 2(0.2)(0.96) (0.2)^2 =0.
Average slope:
0 1 1 2 2 1
f t y f t y
= ⎣^ ⎦=
Correction: y 2 = y 1 + φ 2 h = 0.96 − (0.498511)(0.2) =0.
Exact: y t ( ) = 1/((0.4)^2 + 1) =0.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0
1
Time (sec)
y
Hw12, Problem 1(c)
max error = 0.
Modified-Euler, h=0. Euler, h = 0. Exact
Figure 5. 2 nd^ -order RK with h = 0.2—redline is almost completely hidden by the exact soln.
0 0.05 0.1 0.15 0.2 0.
1
Time (sec)
y
Hw12, Problem 1(c) Modified-Euler, h=0. Euler, h = 0. Exact
Figure 6. Close-up of Figure 5
As seen in Figure 5, the 2 nd^ -order RK with h =0.2 is more accurate than the Euler integration with h = 0.1. The factor of improvement is 6.4. This is even though both methods utilize the exact same number of function calls! Since the modified-Euler has errors of order O(h 2 ), one might anticipate that the modified- Euler would be 0.1/(0.2) 2 = 2.5 times more accurate than Euler integration.
The following Matlab code generates the 2nd^ -order RK solution, the 4th^ -order RK solution, and the ode solution and compares them against the exact:
y_exact_d = 1./(Td.^2 +1);
max_error_d = max(abs(y_exact_d - Yd));
figure(4) plot(Td,Yd,'r-',Ta,Ya,'o--',Td,y_exact_d,'k-','linewidth',2); xlabel('Time (sec)'); ylabel('y'); axis([0 2 0 1]) legend('4th-order RK, h=0.4','Euler, h = 0.1','Exact') gtext(['max error = ',num2str(max_error_d)]); title('Hw12, Problem 1(d)')
% part (d), Using ode45: % y0 = 1; [Tm,Ym] = ode45(f,[0 2],1);
y_exact_m = 1./(Tm.^2 +1);
max_error_d = max(abs(y_exact_m - Ym));
figure(5) plot(Tm,Ym,'r-',Ta,Ya,'o--',Tm,y_exact_m,'k-','linewidth',2); xlabel('Time (sec)'); ylabel('y'); axis([0 2 0 1]) legend('ode45','Euler, h = 0.1','Exact') gtext(['max error = ',num2str(max_error_d)]); title('Hw12, Problem 1(d)')
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0
1
Time (sec)
y
Hw12, Problem 1(d)
max error = 2.1715e-
ode Euler, h = 0. Exact
Figure 7.
Problem 2. (5 points) Find the step response of a spring-mass-damper system having the following governing equation:
2 2 ( )
d x d x m c k x F t d t d t
dx x dt
where F(t) is a unit step function given by
; t
; t F(t)
Note that the problem is second-order, since the highest derivative in x is 2. Before proceeding, we need to write an equivalent set of 2 first-order equations. Define y 1 (t) = x(t) and y 2 (t) = dx/dt. Then the system of equations is:
y f(t,y ) y (F(t) ky cy )/m
y y ⇒ = = − −
2 1 2
1 2
(a) Find the step response from t = 0 to t = 50 seconds using the Matlab function ode45.m. Use the parameter values m = 1kg, k = 1 N/m, c = 0.2 Ns/m. Your call will have the following form:
[T,X] = ode45(func,[0 50],[0 0]);
where func is an anonymous function that returns the function f(t,y). Note that the outputs from ode45 are an Nx1 vector of times, T, and an Nx2 matrix of response, X. The first column of X is the displacement and the second is the velocity. Compare the displacement calculated with ode45 with the exact solution given by:
x( t) = − e −^ nt sin ω dt cos ω dt
where ς = c / ( 2 k/m ), ω n = k/m , ω d = ω n 1 − ς^2
(b) Compare the answer from ode45 with one obtained using first-order Euler integration. Assuming that the number of time steps used by ode45 is equal to N_RK = length(T)-1, use N_E = 4*N_RK steps in the Euler integration scheme. Note that this balances the number of calls to the derivative subroutine between the two methods.
Solution
The next plot compares the exact solution with the ode45 solution. There are actually two curves plotted, but they lie almost on top of one another. The maximum error in the displacement is only 6.02x10 -4. The second plot shows a close-up around a time of 35 seconds.
k
m
c
x(t)
F
0 5 10 15 20 25 30 35 40 45 50
0
1
Time (sec)
y
Hw12, Problem 2(b)
1st-order Euler Exact
It is seen that the RK4 solution is much, much, more accurate that the Euler-integration, even though both methods make the same number of function calls. !!
Matlab Code:
% m-file solves problem 2 of hw % % part (a), ode45: %
% Assign parameters m=1; k=1; c=0.2; F=1;
func = @(t,y) [y(2); (F-ky(1)-cy(2))/m];
[T,X] = ode45(func,[0 50],[0 0]); y = X(:,1);
% calculate exact solution
wn=sqrt(k/m); zeta = c/2/sqrt(km); wd = wnsqrt(1-zeta^2); N = length(T); y_exact = zeros(1,N);
for i = 1:L; t = T(i); term0 = (cos(wdt) + zetawnsin(wdt)/wd); y_exact(i) = 1 - exp(-zetawnt)*term0; end y_exact = y_exact(:); % makes it a column vector max_error_a = max(abs(y_exact - y));
figure(12) plot(T,y,'--',T,y_exact,'r-','linewidth',2); xlabel('Time (sec)') ylabel('y') legend('ode45','Exact') gtext(['max error = ',num2str(max_error_a)]); title('Hw12, Problem 2(a)')
% part (b), Euler with h = 1/4th of RK: % h = 50/N/4; Y_euler = [0 0]; T_euler = 0; tn = 0; yn = [0 0]'; while tn <=50-h;
ynp1 = yn + func(tn,yn)*h; tn = tn + h; yn = ynp1; T_euler = [T_euler tn]; Y_euler = [Y_euler ; yn']; end x_euler = Y_euler(:,1);
figure(13) plot(T_euler,x_euler,'--',T,y_exact,'r-','linewidth',2); xlabel('Time (sec)') ylabel('y') legend('1st-order Euler', 'Exact') title('Hw12, Problem 2(b)')
