LAB02 [Part 1].pdf

42 MATLAB sessions: Laboratory 2 Laboratory 2 Numerical Solutions by Euler and Improved Euler Methods (scalar equations) (x2.4 and x2.5 in the Edwards/Penney text) In this session we look at basic numerical methods to help us understand the fundamentals of numerical approximations. Our objective is as follows. 1. Implement Euler?s method as well as an improved version to numerically solve an IVP. 2. Compare the accuracy and e?ciency of the methods with methods readily available in MATLAB. 3. Apply the methods to speciflc problems and investigate potential pitfalls of the methods. Euler?s Method To derive Euler?s method start from y(t0) = y0 and consider a Taylor expansion at t1 = t0 +h: y(t1) = y(t0)+y0(t0)(t1 ¡t0)+::: = y0 +hf(t0;y(t0))+::: = y0 +hf(t0;y0)+::: For small enough h we get an approximation y1 for y(t1) by suppressing the :::, namely y1 = y0 +hf(t0;y0): (L2.1) Note that the IVP y0 = f(t;y), y(t0) = y0 can be reformulated as y(t1) = y(t0)+ Z t1 t0 f(s;y(s))ds: (L2.2) Then (L2.1) amounts to approximating the integral in (L2.2) using a left point rule. The iteration (L2.1) is repeated to obtain y2 ? y(t2);::: such that yn+1 = yn +hf(tn;yn) tn+1 = tn +h Geometrically, the approximation made is equivalent to replacing the solution curve by the tangent line at (t0;y0). From the flgure we have f(t0;y0) = f(t0;y(t0)) = y0(t0) = tan = y1 ¡y0h ; from which (L2.1) follows. .... ...... ... s s y0 y1 y(t1) t0 t1 h ¡ ¡ ¡ ¡ ¡ ¡ As an example consider the IVP y0 = 2y = f(t;y) with y(0) = 3: Note that here f does not explicitly depend on t (the ODE is called autonomous), but does implicitly through y = y(t). To apply Euler?s method we start with the initial condition and select a step size h. Since we are constructing arrays t and y without dimensionalizing them flrst it is best to clear these names in case they have been used already in the same MATLAB work session. MATLAB sessions: Laboratory 2 43 >> clear t y % no comma between t and y! type help clear for more info >> y(1)=3; t(1)=0; h=0.1; Since f is simple enough we may use the inline syntax: >> f=inline(?2*y?,?t?,?y?) f = Inline function: f(t,y) = 2*y Note that the initialization y(1)=3 should not be interpreted as \the value of y at 1 is 3", but rather \the flrst value in the array y is 3". In other words the 1 in y(1) is an index, not a time value! Unfortunately, MATLAB indices in arrays must be positive (a legacy from Fortran...). The successive approximations at increasing values of t are then obtained as follows: >> y(2)=y(1)+h*f(t(1),y(1)), t(2)=t(1)+h, y = 3.0000 3.6000 t = 0 0.1000 >> y(3)=y(2)+h*f(t(2),y(2)), t(3)=t(2)+h, y = 3.0000 3.6000 4.3200 t = 0 0.1000 0.2000 >> y(4)=y(3)+h*f(t(3),y(3)), t(4)=t(3)+h, y = 3.0000 3.6000 4.3200 5.1840 t = 0 0.1000 0.2000 0.3000 >> y(5)=y(4)+h*f(t(4),y(4)), t(5)=t(4)+h, y = Columns 1 through 4 3.0000 3.6000 4.3200 5.1840 Column 5 6.2208 t = Columns 1 through 4 0 0.1000 0.2000 0.3000 Column 5 0.4000 >> y(6)=y(5)+h*f(t(5),y(5)), t(6)=t(5)+h, y = Columns 1 through 4 3.0000 3.6000 4.3200 5.1840 Columns 5 through 6 6.2208 7.4650 t = Columns 1 through 4 0 0.1000 0.2000 0.3000 Columns 5 through 6 0.4000 0.5000 The arrays y and t are now 1 £ 6 row arrays. The dimension increases as new values of y and t are computed. Each time a new value is computed the whole array is output. To avoid this the output in each command can be suppressed (with a ;) and the list of computed y values vs t values can be output at the end only. To do this we simply concatenate the column versions of t and y into a 6£2 array: