Coursehero >> California >> UCSD >> MAE 107

Course Hero has millions of student submitted documents similar to the one below including study guides, homework solutions, papers, exam answer keys and textbook solutions.

#9 Homework PROBLEM 1. (15 points) Consider the initial value problem: y = 2t 2 - y , y (0) = -1 (1) The exact solution is y (t ) = -5exp(-t ) + 2t 2 - 4t + 4 (2) The objective here is to test various numerical methods for computing the solutions of ordinary differential equations in terms of their efficiency in obtaining an accurate solution. Implement each of the following numerical methods in Matlab to solve (1) with a uniform step size h and obtain y(0.5) accurate to 4 decimal places. a. Explicit Euler b. Modified Euler predictorcorrector c. Fourthorder RungeKutta Which of these methods is most efficient (i.e., requires the fewest evaluations of the righthand side of (1)) ? Solution: Transforming the procedures of examples 20.1, 20.2 and 20.3 into a MATLAB code, the solution of the problem is presented as following: % % % % % % % % % % % MAE 107 S2007, homework 9, problem 1 Exact solution to double precision accuracy y(0.5) = -0.532653298563167 Exact solution to 4 decimal digit accuracy y(0.5) = -0.5326 Euler Method y(0.5) = -0.5326 function evaluations = 512 Modified Euler Method y(0.5) = -0.5326 function evaluations = 64 4th-Order Runge-Kutta Method y(0.5) = -0.5326 function evaluations = 8 function [] = hw9_prob1() clc;clear; % Solution for hw9, problem #1 % integration interval [a b] a = 0; b = 0.5; % initial condition y0 = -1; % number of decimal digit accuracy required Nd = 4; %********************************************* % the exact solution evaluated at t = b format long; yexact = -5*exp(-b) + 2*b^2 - 4*b + 4; disp([char(10) 'Exact solution to double precision accuracy']) disp(['y(' num2str(b) ') = ' num2str(yexact,15)]) yexact = fixdec(yexact,Nd); disp(['Exact solution to ' num2str(Nd) ' decimal digit accuracy']) disp(['y(' num2str(b) ') = ' num2str(yexact,Nd)]) %********************************************* % the explicit Euler method h = 2*(b - a); err = 1; while err > 0 h = h/2; [t, y] = odeEuler(@hw9_prob1_ODE,a,b,h,y0); tN = t(length(t)); yN = fixdec(y(length(t)),Nd); err = abs(yN - yexact); end disp([char(10) 'Euler Method']) disp(['y(' num2str(tN) ') = ' num2str(yN,Nd)]) N = (b - a)/h; % number of steps Nf = N; % number of function evaluations disp(['function evaluations = ', num2str(Nf)]) %********************************************** % the modified Euler method h = 2*(b - a); err = 1; while err > 0 h = h/2; [t, y] = odeModEuler(@hw9_prob1_ODE,a,b,h,y0); tN = t(length(t)); yN = fixdec(y(length(t)),Nd); err = abs(yN - yexact); end disp([char(10) 'Modified Euler Method']) disp(['y(' num2str(tN) ') = ' num2str(yN,Nd)]) N = (b - a)/h; % number of steps Nf = 2*N; % number of function evaluations disp(['function evaluations = ', num2str(Nf)]) %************************************************ % the 4th-order Runge-Kutta method h = 2*(b - a); err = 1; while err > 0 h = h/2; end disp([char(10) '4th-Order Runge-Kutta Method']) disp(['y(' num2str(tN) ') = ' num2str(yN,Nd)]) N = (b - a)/h; % number of steps Nf = 4*N; % number of function evaluations disp(['function evaluations = ', num2str(Nf)]) end function y = fixdec(x,n) % fix to specified number of decimals input: % x = array to be fixed n = integer number of decimals (scalar) % output: % y = x fixed to n decimals f = 10^n; y = fix(x*f)/f; end function dydt = hw9_prob1_ODE(t,y) dydt = 2*t.^2 - y; end function [x, y] = odeEuler(dydx,a,b,h,y0) % odeEuler solves a first order initial value ODE using Euler's explicit % method. Input variables: dydx Name of a function file that calculates % dy/dx. a The first value of x. b The last value of x. h % Step size. y0 The value of the solution y at the first point (initial % value. Output variable: x A vector with the x coordinate the of % solution points. y A vector with the y coordinate of the solution % points. x(1) = a; y(1) = y0; N = (b-a)/h; for i = 1:N x(i+1) = x(i) + h; y(i+1) = y(i) + dydx(x(i),y(i))*h; end end function [x, y] = odeModEuler(dydx,a,b,h,y0) % odeModEuler solves a first order initial value ODE using the modified % Euler method. Input variables: dydx Name of a function file that % calculates dy/dx. a The first value of x. b The last value of % x. h Step size. y0 The value of the solution y at the first % point (initial value. Output variable: x A vector with the x % coordinate of the solution points. y A vector with the y coordinate % of the solution points. x(1) = a; y(1) = y0; N = (b-a)/h; for i = 1:N [t, y] = odeRK4(@hw9_prob1_ODE,a,b,h,y0); tN = t(length(t)); yN = fixdec(y(length(t)),Nd); err = abs(yN - yexact); x(i+1) = x(i) + h; slope0 = dydx(x(i),y(i)); y1 = y(i) + slope0*h; slope1 = dydx(x(i+1),y1); y(i+1) = y(i) + (slope0 + slope1)*h/2; end end function [x, y] = odeRK4(dydx,a,b,h,y0) % odeRK4 solves a first order initial value ODE using Ronge-Kutta fourth % order method. Input variables: dydx Name of a function file that % calculates dy/dx. a The first value of x. b The last value of % x. h Step size. y0 The value of the solution y at the first % point (initial value). Output variable: x A vector with the x % coordinate of the solution points. y A vector with the y coordinate % of the solution points. x(1) = a; y(1) = y0; n = (b-a)/h; for i = 1:n x(i+1) = x(i) + h; K1 = dydx(x(i),y(i)); xhalf = x(i) + h/2; yK1 = y(i) + K1*h/2; K2 = dydx(xhalf,yK1); yK2 = y(i) + K2*h/2; K3 = dydx(xhalf,yK2); yK3 = y(i) + K3*h; K4 = dydx(x(i+1),yK3); y(i+1) = y(i) + (K1 + 2*K2 + 2*K3 + K4)*h/6; end end Output: Exact solution to double precision accuracy y(0.5) = 0.532653298563167 Exact solution to 4 decimal digit accuracy y(0.5) = 0.5326 Euler Method y(0.5) = 0.5326 function evaluations = 512 Modified Euler Method y(0.5) = 0.5326 function evaluations = 64 4thOrder RungeKutta Method y(0.5) = 0.5326 function evaluations = 8 >> It is clear that the RungeKutta method is more efficient in terms of minimum number of RHS function evaluations. PROBLEM 2. (15 points) A flexible cable of uniform density is suspended between two points. The shape of the cable, y(x), is governed by the differential equation: d y dy =C 1 2 dx dx 2 where C is a constant equal to the ratio of the weight per unit length of the cable to the magnitude of the horizontal component of tension in the cable at its lowest point. The cable hangs between two points specified by y(0) = 15m and y(20) = 10m, and C = 0.041 m1. Use Matlab's builtin functions ode45 and fzero to devise a shooting method to determine the shape of the cable between x = 0 and x = 20m. Solution: % MAE 107 Spring 2007 % main function to compute shape of flexible cable with appropriate % boundary conditions satisfied function hw9_prob2 clear; close all; clc % setup constants and parameters C = 0.041; xl = 0; xr = 20; yl = 15; yr = 10; % subfunction describing shape of the cable function [dy] = cable_func(x, y) dy = [ y(2); C * sqrt(1 + y(2).^2) ]; end % subfunction to compute residual of boundary value of y at x = xr function resid = newton_func(y_dash_init) y_init = [yl; y_dash_init]; [X, Y] = ode45(@cable_func, [xl xr], y_init); resid = Y(end,1) yr; end % compute that value of y_dash at x = xl for which residual vanishes y_dash_init = fzero(@newton_func, 0); 2 fprintf('The value of dy/dx at x = %8.6f, is %8.6f\n\n', xl, y_dash_init) % plot the shape of the cable y_init = [yl; y_dash_init]; [X, Y] = ode45(@cable_func, [xl xr], y_init); plot(X, Y(:,1)); title('Shape of flexible cable'); end Output: The value of dy/dx at x = 0.000000, is 0.697719

Find millions of documents here - Study Guides, Homework Solutions, Papers, Exam Answer Keys and more. Course Hero has millions of course related materials that will enable you to learn better, faster and get an A in all your courses.
Below is a small sample set of documents: