Ma 449: Numerical Applied Mathematics Model Solutions to Homework Assignment 9 Prof. Wickerhauser Exercise 4 of Section 7.3, p.390. First we note that the integral of t^k from 0 to 4 is 4^(k+1)/(k+1) and takes the following values: k integral of t^k from 0 to 4 - --------------------------- 0 4 1 8 2 64/3 3 64 4 1024/5 Next, set the values of these integrals equal to a linear combination of function values 0^k, 1^k, 2^k, 3^k, and 4^k, with unknown weights w0,w1,w2,w3,w4, and solve for the w's: MATLAB CODE: t=[0 1 2 3 4]; a=[ones(size(t)); t; t.*t; t.*t.*t; t.*t.*t.*t]; b=[4; 8; 64/3; 64; 1024/5]; format long; w=a\b % W = % .311111111111114 % 1.422222222222215 % .533333333333342 % 1.422222222222218 % .311111111111112 (45/2)*w % ANS = % 7.000000000000060 % 31.999999999999833 % 12.000000000000185 % 31.999999999999897 % 7.000000000000024 Except for the round-off error, we recognize that W = (2/45)*[ 7 32 12 32 7 ] Algorithm 3 of Section 7.3, p.392. We do this by Romberg integration: C PROGRAM: /* Section 7.3, Algorithm 3, done with Romberg integration */ #include #include #include #define PI 3.14159265358979323 #define JMAX 15 /* Recursion depth limit. */ /* Function to integrate to get Phi: */ double f(double t) { return exp(-t*t); } /* Recursive composite trapezoid rule increment on 2^J segments of [a,b]: */ double rcti(double a, double b, int j) { double h, sum; int n, k; assert(a=0); assert(j<30); n = 1< #include #include int depth; /* used to track recursion depth */ double f(double x){ return 1/sqrt(x); } double simpson(double a, double b) { return (b-a)*(f(a)+4*f((a+b)/2)+f(b))/6; } double rsrule(double a, double b, double sum, double tol, int descend) { double sr1, sr2, c, err; if(depth>descend) depth=descend; /* track recursion depth */ if(!descend) return sum; /* recursion limit is reached */ c = (a+b)/2; /* else descend if necessary */ sr1 = simpson(a,c); sr2 = simpson(c,b); err = fabs(sr1+sr2 - sum); /* test the current error */ if(err #include #include int depth; /* used to track recursion depth */ double f(double x){ return 1/sqrt(x); } double trapezoid(double a, double b) { return (b-a)*(f(a)+f(b))/2; } double rtrule(double a, double b, double sum, double tol, int descend) { double tr1, tr2, c, err; if(depth>descend) depth=descend; /* track recursion depth */ if(!descend) return sum; /* recursion limit is reached */ c = (a+b)/2; /* else descend if necessary */ tr1 = trapezoid(a,c); tr2 = trapezoid(c,b); err = fabs(tr1+tr2 - sum); /* test the current error */ if(err