s09 - Ma 449: Numerical Applied Mathematics Model Solutions...

Info iconThis preview shows pages 1–4. Sign up to view the full content.

View Full Document Right Arrow Icon

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: 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---------------------------- 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(t);t;t.*t;t.*t.*t;t.*t.*t.*t]; <>b=[4;8;64/3;64;1024/5]; <>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 <assert.h> #include <math.h> #include <stdio.h> #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<b); assert(j>=0); assert(j<30); n = 1<<j; h = (b-a)/n; if(j==0) sum = (f(a)+f(b))/2; else { sum = 0; for(k=1; k<n; k+=2) sum += f(a+k*h); } return h*sum; } main(void) { double r[JMAX][JMAX], phi; int i, j, k; double a=0.0, b[8]={0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0}; double err, tol=5.0e-14; for(i=0; i<8; i++) { r[0][0] = rcti(a,b[i],0); /* Initial trapezoid rule value */ for(j=1; j<JMAX; j++) { r[j][0]=0.5*r[j-1][0]+rcti(a,b[i],j); /* recursive trapezoid rule */ for(k=1; k<=j; k++) /* Romberg integration */ r[j][k] = r[j][k-1] + (r[j][k-1] - r[j-1][k-1])/((1<<(2*k))-1.0); err = fabs( r[j][j] - r[j-1][j-1] ); phi = 0.5 + r[j][j]/sqrt(2*PI); /* adjustment, to get Phi */ printf("phi[%3.1f]=%20.14f err=%8.2e, J=%d\n", b[i], phi, err, j ); if(err<tol) break; /* escape j-loop, go to next value b[i] */ } } return 0; } OUTPUT: phi[0.5]= 0.69147333953143 err=9.33e-03, J=1 phi[0.5]= 0.69146244243255 err=2.73e-05, J=2 phi[0.5]= 0.69146246128594 err=4.73e-08, J=3 phi[0.5]= 0.69146246127401 err=2.99e-11, J=4 phi[0.5]= 0.69146246127401 err=6.05e-15, J=5 phi[1.0]= 0.84152905199630 err=5.28e-02, J=1 phi[1.0]= 0.84134391691401 err=4.64e-04, J=2 phi[1.0]= 0.84134474699460 err=2.08e-06, J=3 phi[1.0]= 0.84134474606877 err=2.32e-09, J=4 phi[1.0]= 0.84134474606854 err=5.82e-13, J=5 phi[1.0]= 0.84134474606854 err=3.33e-16, J=6 phi[1.5]= 0.93325240117164 err=9.25e-02, J=1 phi[1.5]= 0.93320473311640 err=1.19e-04, J=2 phi[1.5]= phi[1....
View Full Document

Page1 / 6

s09 - Ma 449: Numerical Applied Mathematics Model Solutions...

This preview shows document pages 1 - 4. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online