Hw5_solns - logx=log(Panelvals logTrap=log(traperrs logSimp=log(simperrs trapFit=polyfit(logx,logTrap,1 simpFit=polyfit(logx,logSimp,1

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

View Full Document Right Arrow Icon
############ # Solutions for HW5 # by J.R. Gladden ############# from pylab import * ###### #Prob 1 ####### def trap(f,xmin,xmax,numPanels): Isum=0. h=(xmax-xmin)/float(numPanels) for i in range(numPanels): x=xmin+i*h Isum+=(f(x+h)+f(x))*h/2. return Isum def simp13(f,xmin,xmax,numPanels): Isum=0. h=(xmax-xmin)/float(numPanels) for i in range(0,numPanels,2): x=xmin+i*h Isum+=(f(x)+4*f(x+h)+f(x+2*h))*h/3 return Isum def f(x): return sin(x) Itrue=2.0 Panelvals=range(10,500,20) traperrs=zeros(len(Panelvals)) simperrs=zeros(len(Panelvals)) for i in range(len(Panelvals)): Itrap=trap(f,0,np.pi,Panelvals[i]) traperrs[i]=abs(Itrap-Itrue) Isimp=simp13(f,0,np.pi,Panelvals[i]) simperrs[i]=abs(Isimp-Itrue) loglog(Panelvals,traperrs,'b-o',label='Trapezoid') loglog(Panelvals,simperrs,'r-s',label='Simpson') xlabel('Number of Panels',fontsize=18) ylabel('Error',fontsize=18) legend() ############## #Prob 2 #############
Background image of page 1

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

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

Unformatted text preview: logx=log(Panelvals) logTrap=log(traperrs) logSimp=log(simperrs) trapFit=polyfit(logx,logTrap,1) simpFit=polyfit(logx,logSimp,1) print "For Trapeziod method, error scales like N to the: ",trapFit[0] print "For Simpson method, error scales like N to the: ",simpFit[0] ###### #Prob 3 ###### def richardson_trap(g,xmin,xmax,numPanels): ''' Performs a Richardson Extrapolation step using h2=h1/2 and error ~ h^2 ''' p=2 return (2.0**p*trap(g,xmin,xmax,numPanels*2) - trap(g,xmin,xmax,numPanels))/ (2.0**p-1) def g(x): return x*sin(x) Itrue=np.pi N=10 Itrap=trap(g,0,np.pi,N) Irich=richardson_trap(g,0,np.pi,N) print '='*20 print 'Itrue = %2.10f' % Itrue print 'Itrap = %2.10f \t error = %2.4e ppm' % (Itrap,abs(Itrap-Itrue)/Itrue*10**6) print 'Irich = %2.10f \t error = %2.4e ppm' % (Irich,abs(Irich-Itrue)/Itrue*10**6)...
View Full Document

This note was uploaded on 10/05/2010 for the course PHYS phy503 taught by Professor Gladden during the Spring '09 term at Ole Miss.

Page1 / 2

Hw5_solns - logx=log(Panelvals logTrap=log(traperrs logSimp=log(simperrs trapFit=polyfit(logx,logTrap,1 simpFit=polyfit(logx,logSimp,1

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

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