Ma 449: Numerical Applied Mathematics Model Solutions to Homework Assignment 7 Prof. Wickerhauser Exercise 4 of Section 6.1, p.335. (a) Formula 10 MATLAB code and its output: num = [1 -8 0 8 -1]; dx=[-2 -1 0 1 2]; x1=2.3*ones(size(dx)) +0.1*dx, f1=exp(x1) % X1 = 2.1000 2.2000 2.3000 2.4000 2.5000 % F1 = 8.1662 9.0250 9.9742 11.0232 12.1825 x01=2.3*ones(size(dx)) +0.01*dx, f01=exp(x01) % X01 = 2.2800 2.2900 2.3000 2.3100 2.3200 % F01 = 9.7767 9.8749 9.9742 10.0744 10.1757 h=0.1; den=12*h; df1=(f1*num')/den; format long; df1, format short e; err=df1-exp(2.3) % DF1 = 9.974149167936689 ERR = -3.3287E-05 h=0.01; den=12*h; df01=(f01*num')/den; format long; df01, format short e; err=df01-exp(2.3) % DF01 = 9.974182451489710 ERR = -3.3250E-09 (b) Formula 29 MATLAB code and its output: x0=2.3; h=0.1; d2=(exp(x0+2*h)-exp(x0-2*h))/(4*h) % D2 = 10.0408 d1=(exp(x0+h)-exp(x0-h))/(2*h) % D1 = 9.9908 format long; fe=(4*d1-d2)/3, format short e; err=fe-exp(x0) % FE = 9.974149167936696 ERR = -3.3287E-05 h=0.01; d2=(exp(x0+2*h)-exp(x0-2*h))/(4*h) % D2 = 9.9748 d1=(exp(x0+h)-exp(x0-h))/(2*h) % D1 = 9.9743 format long; fe=(4*d1-d2)/3, format short e; err=fe-exp(x0) % FE = 9.974182451489680 ERR = -3.3250E-09 (c) Using |f'''''(c)| < 12.2, we get |E_d| < h^4 * 12.2 / 30, which is evaluated by MATLAB as follows: format short e; h=0.1; h4=h*h*h*h; f5=12.2; ed=h4*f5/30 % ED = 4.0667E-05 h=0.01; h4=h*h*h*h; f5=12.2; ed=h4*f5/30 % ED = 4.0667E-09 Exercise 16 of Section 6.1, p.338. (a) Formula 22 in MATLAB, and its output: num = [1 -8 0 8 -1]; dh=[-2 -1 0 1 2]; x0=3; h=0.1; x=x0*ones(size(num))+h*dh; f=[1.02962 1.06471 1.09861 1.13140 1.16315]; format long; df=(f*num')/(12*h), format short e; err= df-1/3 % DF = .333324999999999 ERR = -8.3333E-06 h=0.001; x=x0*ones(size(num))+h*dh; f=[1.09795 1.09828 1.09861 1.09895 1.09928]; format long; df=(f*num')/(12*h), format short e; err= df-1/3 % DF = .335833333333443 ERR = .0025 (b) Formula 24 can be estimated as follows: E = E_r + E_d, so |E| <= |E_r| + |E_d|. Now for f(x)=ln(x), we have f''''(x) = 24/x^5, which is no bigger than 24/(x0-2*h)^5 on the interval x-2h #include #define MAX1 15 /* Set which part of the problem we will solve: */ #define X0 X0a; #define f fa; #define X0e 0.0001; double fe(double x){ return pow(x,pow(x,x));} #define X0d (1-sqrt(5.))/2; double fd(double x){ return sin(pow(x,3.)-7*pow(x,2.)+6*x+8);} #define X0c 1/sqrt(2.); double fc(double x){ return sin(cos( 1/x));} #define X0b (1+sqrt(5.))/3; double fb(double x){ return tan(cos((sqrt(5.)+sin(x))/(1+x*x)));} #define X0a (1/sqrt(3.)) double fa(double x){ return 60*pow(x,45.)-32*pow(x,33.)+233*pow(x,5.)-47*x*x-77.; } void message(const char *text) { printf( "%s.\n", text); } /* O(h^4) derivative formula: */ double df(double x, double h) { return ( -f(x+2*h)+8*f(x+h)-8*f(x-h)+f(x-2*h) ) / ( 12*h ); } int main(void) { double h, x, tolerance, eps, D[MAX1], H[MAX1], R[MAX1], E[MAX1]; int k; eps = 1.0e-15; /* Add this to avoid relative error problems */ tolerance = 1.0e-13; /* Target accuracy for the derivative */ x=X0; /* Point at which to differentiate */ h=1; /* Initial h */ D[0]=df(x,h); /* First approximation to the derivative */ E[0]=R[0]=0; /* Initialize error trackers */ for(k=1; 1 ; k++) { h /= 10; /* Decrease h for the next trial */ D[k]=df(x,h); E[k]=fabs(D[k]-D[k-1]); R[k]=2*E[k] / ( fabs(D[k]) + fabs(D[k-1]) + eps ); if ( R[k]1 && E[k]>E[k-1]) { message("Total error has begun to increase"); break; } if( k==MAX1-1 ) { message("Maximum number of iterations reached"); break; } } printf("Iteration %d: h=%6.1e, f'(x0)= %20.15e, E=%6.1e, R=%6.1e\n", k, h, D[k], E[k], R[k] ); return 0; } MATLAB CODE: % df4.m: % Fourth-order centered difference approximation to the first derivative function df = df4(f, x, h) df=(-f(x+2*h)+8*f(x+h)-8*f(x-h)+f(x-2*h))./(12*h); % for one x, many h end % Main computation: h=0.1.^(0:14); format short g; h % h = 1 0.1 0.01 ... 1e-14 % Part (a) function and evaluation point: x = 1/sqrt(3.); f = @(x) 60*x.^45 - 32*x.^33 + 233*x.^5 - 47*x.*x - 77; % common computation of derivatives and errors: D = df4(f, x, h); E = abs(diff(D)); R = 2*E./(abs(D(1:14))+abs(D(2:15))); format long; [h(2:15)' D(2:15)' E' R'] % ans = % 1.00000000000000e-01 7.50857250132671e+01 1.59137417504566e+19 2.00000000000000e+00 % 1.00000000000000e-02 7.51734854395347e+01 8.77604262675504e-02 1.16812042340815e-03 % 1.00000000000000e-03 7.51734946942512e+01 9.25471655932597e-06 1.23111439299821e-07 % 1.00000000000000e-04 7.51734946952235e+01 9.72264047049975e-10 1.29336018100279e-11 % 1.00000000000000e-05 7.51734946944064e+01 8.17124146124115e-10 1.08698438118572e-11 % 1.00000000000000e-06 7.51734947073146e+01 1.29081882960236e-08 1.71711962406290e-10 % 1.00000000000000e-07 7.51734946173125e+01 9.00020893368492e-08 1.19725828553198e-09 % 1.00000000000000e-08 7.51734957541809e+01 1.13686837721616e-06 1.51232608568627e-08 % 1.00000000000000e-09 7.51734875829394e+01 8.17124144703030e-06 1.08698442305398e-07 % 1.00000000000000e-10 7.51733727118638e+01 1.14871075624023e-04 1.52808080459787e-06 % 1.00000000000000e-11 7.51728161200541e+01 5.56591809697693e-04 7.40413593483164e-06 % 1.00000000000000e-12 7.51801583949903e+01 7.34227493619244e-03 9.76671723306381e-05 % 1.00000000000000e-13 7.51991062012772e+01 1.89478062869313e-02 2.52000252000189e-04 % 1.00000000000000e-14 7.55543775691573e+01 3.55271367880036e-01 4.71327572662982e-03 % Part (b) function and evaluation point: x = (1+sqrt(5.))/3; f = @(x) tan(cos((sqrt(5.)+sin(x))./(1+x.*x))); % common computation of derivatives and errors: D = df4(f, x, h); E = abs(diff(D)); R = 2*E./(abs(D(1:14))+abs(D(2:15))); format long; [h(2:15)' D(2:15)' E' R'] % ans = % 1.00000000000000e-01 1.22871864901577e+00 2.06670436685430e-02 1.69626526190058e-02 % 1.00000000000000e-02 1.22859743490954e+00 1.21214106231893e-04 9.86556894530764e-05 % 1.00000000000000e-03 1.22859742337716e+00 1.15323770533138e-08 9.38661988032549e-09 % 1.00000000000000e-04 1.22859742337900e+00 1.83653092733493e-12 1.49481912658201e-12 % 1.00000000000000e-05 1.22859742340368e+00 2.46793696589975e-11 2.00874340032996e-11 % 1.00000000000000e-06 1.22859742311456e+00 2.89120505314600e-10 2.35325664730474e-10 % 1.00000000000000e-07 1.22859742402818e+00 9.13620734621645e-10 7.43629049754858e-10 % 1.00000000000000e-08 1.22859740529317e+00 1.87350137625941e-08 1.52491072657584e-08 % 1.00000000000000e-09 1.22859764029037e+00 2.34997207249066e-07 1.91272734064347e-07 % 1.00000000000000e-10 1.22859676136381e+00 8.78926561309612e-07 7.15390333559215e-07 % 1.00000000000000e-11 1.22859754777179e+00 7.86407975628123e-07 6.40086111793713e-07 % 1.00000000000000e-12 1.22860286759045e+00 5.31981865958819e-06 4.32998352623504e-06 % 1.00000000000000e-13 1.22950261083332e+00 8.99743242873452e-04 7.32062355151986e-04 % 1.00000000000000e-14 1.21638810135494e+00 1.31145094783849e-02 1.07237084739995e-02 % Part (c) function and evaluation point: x = 1/sqrt(2.); f = @(x) sin(cos( 1./x)); % common computation of derivatives and errors: D = df4(f, x, h); E = abs(diff(D)); R = 2*E./(abs(D(1:14))+abs(D(2:15))); format long; [h(2:15)' D(2:15)' E' R'] % ans = % 1.00000000000000e-01 1.96076386074356e+00 9.32301345020599e-01 6.23774333292646e-01 % 1.00000000000000e-02 1.95156036922614e+00 9.20349151741995e-03 4.70487156811705e-03 % 1.00000000000000e-03 1.95155960901272e+00 7.60213421768796e-07 3.89541405853383e-07 % 1.00000000000000e-04 1.95155960893668e+00 7.60409513134164e-11 3.89641961043488e-11 % 1.00000000000000e-05 1.95155960892151e+00 1.51731960329471e-11 7.77490780375051e-12 % 1.00000000000000e-06 1.95155960902351e+00 1.02001962432041e-10 5.22668956474996e-11 % 1.00000000000000e-07 1.95155960761954e+00 1.40397005132797e-09 7.19409258800701e-10 % 1.00000000000000e-08 1.95155962108099e+00 1.34614541735800e-08 6.89779296240549e-09 % 1.00000000000000e-09 1.95155951815407e+00 1.02926926093261e-07 5.27408579761840e-08 % 1.00000000000000e-10 1.95155977026721e+00 2.52113145027266e-07 1.29185467518335e-07 % 1.00000000000000e-11 1.95156182880574e+00 2.05853852452975e-06 1.05481649611874e-06 % 1.00000000000000e-12 1.95152459007512e+00 3.72387306173394e-05 1.90816838885248e-05 % 1.00000000000000e-13 1.95232718880334e+00 8.02598728218529e-04 4.11182992428626e-04 % 1.00000000000000e-14 1.95167955870564e+00 6.47630097698526e-04 3.31777140555281e-04 % Part (d) function and evaluation point: x = (1-sqrt(5.))/2; f = @(x) sin(x.^3 -7*x.^2+6*x+8); % common computation of derivatives and errors: D = df4(f, x, h); E = abs(diff(D)); R = 2*E./(abs(D(1:14))+abs(D(2:15))); format long; [h(2:15)' D(2:15)' E' R'] % ans = % 1.00000000000000e-01 3.67331992907900e+00 3.97719279426030e+00 2.00000000000000e+00 % 1.00000000000000e-02 2.96568581351550e+00 7.07634115563501e-01 2.13174726156197e-01 % 1.00000000000000e-03 2.96551484641162e+00 1.70967103879871e-04 5.76500825659039e-05 % 1.00000000000000e-04 2.96551482918591e+00 1.72257115238494e-08 5.80867487300644e-09 % 1.00000000000000e-05 2.96551482917203e+00 1.38782318970243e-11 4.67987270219338e-12 % 1.00000000000000e-06 2.96551482934874e+00 1.76710646115907e-10 5.95885221588927e-11 % 1.00000000000000e-07 2.96551482731333e+00 2.03540917453893e-09 6.86359466185653e-10 % 1.00000000000000e-08 2.96551485682676e+00 2.95134285899223e-08 9.95221071607277e-09 % 1.00000000000000e-09 2.96551464865994e+00 2.08166816673128e-07 7.01958459254185e-08 % 1.00000000000000e-10 2.96551505574172e+00 4.07081775399831e-07 1.37271872065560e-07 % 1.00000000000000e-11 2.96549359142991e+00 2.14643118101598e-05 7.23799713911924e-06 % 1.00000000000000e-12 2.96547971364210e+00 1.38777878073704e-05 4.67976741540970e-06 % 1.00000000000000e-13 2.96633088462765e+00 8.51170985545657e-04 2.86985220261057e-04 % 1.00000000000000e-14 2.95874436062604e+00 7.58652400160553e-03 2.56081946222801e-03 % Part (e) function and evaluation point: x = 0.0001; f = @(x) x.^(x.^x); % common computation of derivatives and errors: D = df4(f, x, x*h); % scale h to fractions of x to avoid imaginary parts E = abs(diff(D)); R = 2*E./(abs(D(1:14))+abs(D(2:15))); format long; [x*h(2:15)' D(2:15)' E' R'] % show the scaled h values % ans = % 1.00000000000000e-05 1.01520613064893e+00 4.87556088335509e-04 4.80274821829985e-04 % 1.00000000000000e-06 1.01520610590270e+00 2.47462277336297e-08 2.43755699341703e-08 % 1.00000000000000e-07 1.01520610590052e+00 2.18092210957366e-12 2.14825550880308e-12 % 1.00000000000000e-08 1.01520610590065e+00 1.35447209004269e-13 1.33418434165262e-13 % 1.00000000000000e-09 1.01520610588981e+00 1.08419939692794e-11 1.06795988581242e-11 % 1.00000000000000e-10 1.01520610569104e+00 1.98770333525999e-10 1.95793083189979e-10 % 1.00000000000000e-11 1.01520610697853e+00 1.28748989247640e-09 1.26820542591557e-09 % 1.00000000000000e-12 1.01520609749176e+00 9.48676914802604e-09 9.34467309361056e-09 % 1.00000000000000e-13 1.01520557910760e+00 5.18384163594732e-07 5.10619762060222e-07 % 1.00000000000000e-14 1.01520862842621e+00 3.04931861005997e-06 3.00364191576826e-06 % 1.00000000000000e-15 1.01521980926111e+00 1.11808349037013e-05 1.10132765047719e-05 % 1.00000000000000e-16 1.01476805835591e+00 4.51750905202219e-04 4.45077443475092e-04 % 1.00000000000000e-17 1.01971473076788e+00 4.94667241196511e-03 4.86283043376892e-03 % 1.00000000000000e-18 1.00966327312713e+00 1.00514576407513e-02 9.90594913462104e-03 The indented row, which has the smallest R value, contains the best h and the most accurate D value. Exercise 4 of Section 6.2, p.349. NOTE: part b has a typo: Should be "use h=0.1" rather than "use h=0.01". MATLAB Code: (a) num=[1 -2 1]; x0=1; h=0.05; f=[.5817 .5403 .4976] format long; d2f=(f*num')/(h*h), -cos(x0) format short e; err=d2f+cos(x0) % D2F = -.520000000000009 % ANS = -.540302305868140 ERR = .0203 (b) num=[1 -2 1]; x0=1; h=0.1; f=[.6216 .5403 .4536]; format long; d2f=(f*num')/(h*h), -cos(x0) format short e; err=d2f+cos(x0) % D2F = -.539999999999996 % ANS = -.540302305868140 ERR = 3.0231E-04 (c) num=[-1 16 -30 16 -1]; x0=1; h=0.05; f=[.6216 .5817 .5403 .4976 .4536]; format long; d2f=(f*num')/(12*h*h), -cos(x0) format short e; err=d2f+cos(x0) % D2F = -.513333333333358 % ANS = -.540302305868140 ERR = .0270 (c') num=[-1 16 -30 16 -1]; x0=1; h=0.1; f=[.6967 .6216 .5403 .4536 .3624]; format long; d2f=(f*num')/(12*h*h), -cos(x0) format short e; err=d2f+cos(x0) % D2F = -.540833333333323 % ANS = -.540302305868140 ERR = -5.3103E-04 (d) Value b is the most accurate of a,b,c. Value c' is computed to show how the error actually decreases, even for the higher-order formula, when the step size h is increased. Exercise 7 of Section 6.2, p.349. Subtracting: f(x+h)=f(x) + hf'(x) + h^2f''(x)/2 + h^3f'''(x)/6 + h^4f''''(x)/24 +O(h^5) f(x-h)=f(x) - hf'(x) + h^2f''(x)/2 - h^3f'''(x)/6 + h^4f''''(x)/24 +O(h^5) --------------------------------------------------------------------------- f(x+h)-f(x-h) = 2hf'(x) + h^3f'''(x)/3 +O(h^5) Subtracting: f(x+2h)=f(x) + 2hf'(x) + 2h^2f''(x) + 4h^3f'''(x)/3 + 2h^4f''''(x)/3 +O(h^5) f(x-2h)=f(x) - 2hf'(x) + 2h^2f''(x) - 4h^3f'''(x)/3 + 2h^4f''''(x)/3 +O(h^5) --------------------------------------------------------------------------- f(x+2h)-f(x-2h) = 4hf'(x) + 8h^3f'''(x)/3 +O(h^5) Eliminate f' by another subtraction: f(x+2h)-f(x-2h) = 4hf'(x) + 8h^3f'''(x)/3 +O(h^5) 2[f(x+h)-f(x-h)] = 4hf'(x) + 2h^3f'''(x)/3 +O(h^5) --------------------------------------------------------------------------- f(x+2h)-f(x-2h) - 2[f(x+h)-f(x-h)] = 6h^3f'''(x)/3 +O(h^5) Thus, f'''(x) = [f(x+2h) - 2f(x+h) + 2f(x-h) - f(x-2h)]/2h^3 + O(h^2) Algorithm 1 of Section 6.2, p.351. C Code: #include #include #include /* Differentiation formulas based on Newton polynomial interpolation. */ /* Compute the Newton polynomial leading coefficients a[0],...,a[n] using the divided differences derived from {(x[k],f[k]) : k=0,1,...,n}: */ void divided_differences ( double *a, const double *x, const double *f, int n) { int k, j; for(k=0; k<=n; k++) a[k] = f[k]; /* Initialize a[] with f's: */ /* Compute jth divided difference, putting it in a[j]: */ for(j=1; j<=n; j++) for(k=n; k>=j; k--) a[k] = (a[k]-a[k-1])/(x[k]-x[k-j]); return; } /* Use Horner's method: */ double evaluate_Newton_polynomial(const double *a, const double *x, int n, double t) { int k; double y; y = a[n]; for(k=n-1; k>=0; k--) y = a[k] + (t-x[k])*y; return y; } /* Permute x0 <- x1 <- ... <- xn <- x0 */ void cycle( double *z, int n) { double temp; int k; temp = z[0]; for(k=0; k