2071 Lab 02 {@(#) Tue Jan 31 11:35:30 2017 } $Id: summary.txt,v 1.10 2017/03/14 19:05:06 mike Exp $ EXERCISE 1. 1. copy forward_euler.m 2. copy expm_ode.m 3. Euler's explicit method numSteps Stepsize Euler Error Ratio 10 0.2 -3.21475 0.055922 2.03226 20 0.1 -3.24315 0.027517 2.01647 40 0.05 -3.25702 0.013646 2.00829 80 0.025 -3.26388 0.006795 2.00416 160 0.0125 -3.26728 0.003390 2.00208 320 0.00625 -3.26898 0.001693 4. first order. EXERCISE 2. 1. rk2.m: CHECK comments xa = x(1,k)+h/2 ; ya = y(:,k)+h/2*f_ode(x(1,k),y(:,k)) ; y(:,k+1) = y(:,k) + h * f_ode(xa,ya); 2. RK2 numSteps Stepsize RK2 Error Ratio 10 0.2 -3.27490 -4.2255e-03 4.3367e+00 20 0.1 -3.27164 -9.7435e-04 4.1588e+00 40 0.05 -3.27090 -2.3429e-04 4.0771e+00 80 0.025 -3.27073 -5.7464e-05 4.0380e+00 160 0.0125 -3.27068 -1.4231e-05 4.0189e+00 320 0.00625 -3.27067 -3.5409e-06 3. quadratic 4. Euler numSteps Stepsize Error RK2 Error 10 0.2 5.5922e-02 4.2255e-03 20 0.1 2.7517e-02 9.7435e-04 40 0.05 1.3646e-02 2.3429e-04 80 0.025 6.7950e-03 5.7464e-05 160 0.0125 3.3904e-03 1.4231e-05 320 0.00625 1.6935e-03 3.5409e-06 5. Euler at 160 is about the same as rk2 at 10 6. euler(320)=1.69e-3: 8.45e-3, 4.23e-3, 2.12e-3, 1.06e-4, 5.30e-5, 2.65e-5, 1.33e-5, 6.65e-6, 3.32e-6 9 halvings is 163840 steps OR: N/320=1.6935e-3/3.5409e-6: N=153050 7. Check prediction accuracy is close. EXERCISE 3 1. rk3.m all 2. RK3 numSteps Stepsize RK3 Error Ratio 10 0.2 -3.27046 2.11794e-04 8.66704e+00 20 0.1 -3.27065 2.44367e-05 8.32700e+00 40 0.05 -3.27067 2.93463e-06 8.16178e+00 80 0.025 -3.27067 3.59558e-07 8.08045e+00 160 0.0125 -3.27067 4.44973e-08 8.04013e+00 320 0.00625 -3.27067 5.53440e-09 3. cubic 4. numSteps Stepsize RK2 Error RK3 Error 10 0.2 -4.2255e-03 2.11794e-04 20 0.1 -9.7435e-04 2.44367e-05 40 0.05 -2.3429e-04 2.93463e-06 80 0.025 -5.7464e-05 3.59558e-07 160 0.0125 -1.4231e-05 4.44973e-08 320 0.00625 -3.5409e-06 5.53440e-09 5. rk2 takes about 40 steps to achieve the error rk3 gets with 10 steps 6. Crk2=3.54e-06/(320^2)=.3625, n=sqrt(Crk2/5.53e-9)=8096 7. Check prediction accuracy is close. EXERCISE 4 1. fValue = [-8;-9] 2. solve ysystem: ysystem(:,end) = [ -2.7293;-2.5940 ] 3. solve yscalar: same values 4. norm comparison: both are zero using octave and Matlab. EXERCISE 5 1. pendulum_ode.m y1prime=y2, y2prime=-3*sin(y1) 2. x(1)=0, x(251)=6.25, x(501)=12.50, x(751)=18.75, x(1001)=25 3 x n Euler RK3 0.00 1 1 0 1 0 6.25 251 -1.06185 0.97913 -0.75431 1.06326 12.50 501 1.21137 -1.33835 0.11877 -1.64791 18.75 751 -1.82570 0.35516 0.58196 1.33043 25.00 1001 0.90952 2.72745 -0.97403 -0.35971 4. plots: curves close, Euler growing away from RK. BE SURE THEY start from y=1. EXERCISE 6 CHECK dx sizes 1. Euler, [0,20], 40: kinda parabolic down with straight tail, dx=.5 2 Euler [0,20], 30, 20, 15, 12, 10: +/- osc grows, last is const dx(30)=.66667, dx(20)=1, dx(15)=1.3333, dx(12)=1.6667, dx(10)=2 3. Euler, [0,20], 8: explodes, max approx 400, dx=2.5 4. rk3, [0,20], 20, 10, 9, 8: same as 2, first is stable dx(20)=1, dx(10)=2, dx(9)=2.2222, dx(8)=2.5 5. rk3, [0,20], 7: explodes, peak around -650, dx=2.8571 EXERCISE 7 1. ab2.m 2. add comments: % The Adams-Bashforth algorithm starts here % The Runge-Kutta algorithm starts here 3. why is dydxold used? 4. if numSteps=100, f gets called 101 times 5. AB2 numSteps Stepsize AB2(x=2) Error(x=2) Ratio 10 0.2 -3.28014 -9.46936e-03 4.08540e+00 20 0.1 -3.27299 -2.31786e-03 4.05194e+00 40 0.05 -3.27124 -5.72036e-04 4.02808e+00 80 0.025 -3.27081 -1.42012e-04 4.01453e+00 160 0.0125 -3.27071 -3.53745e-05 4.00738e+00 320 0.00625 -3.27068 -8.82734e-06 6. quadratic EXERCISE 8 (Extra Credit: stability regions) 1. m-file to plot unit circle 2. m-file to plot mu=1-zeta=stability region explicit Euler circle centered at -1. 3. slightly smaller circle designates inside in different color. 4. add x-axis, y-axis 5. copy m-file, plot both explicit Euler and AB2 stability (smaller circle squished up against imag axis) 6. Generate ODEs and h values so that (a) both ab2 and explicit Euler stable (b) ab2 stable, not explicit Euler (c) explicit Euler stable, not ab2 [lambda=-2+10i, h=0.06 is one example] EXERCISE 8 (Extra Credit) (OMITTED) 1. heun.m 2. first column of: (final ratio of errors=4.0289) 3. Heun AB2 RK2 numSteps Error Calls Error Calls Error Calls 10 2.0133e-02 20 4.7542e-02 11 4.1526e-03 20 20 5.0656e-03 40 1.2208e-02 21 3.8664e-04 40 40 1.2430e-03 80 3.0518e-03 41 4.0838e-05 80 80 3.0677e-04 160 7.5934e-04 81 4.4799e-06 160 160 7.6144e-05 320 1.8913e-04 161 4.7183e-07 320 r=4.0289 r=4.0150 r approx 4 4. all second order 5. Heun vs AB2: Heun more accurate, ab2 faster measured in calls 6. all 3: rk2 most accurate, fastest in terms of calls