{[ promptMessage ]}

Bookmark it

{[ promptMessage ]}

tutor-notes-2005-sheet-7

tutor-notes-2005-sheet-7 - ”Ofe‘é £v Tim's The...

Info icon This preview shows pages 1–12. Sign up to view the full content.

View Full Document Right Arrow Icon
Image of page 1

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

View Full Document Right Arrow Icon
Image of page 2
Image of page 3

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

View Full Document Right Arrow Icon
Image of page 4
Image of page 5

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

View Full Document Right Arrow Icon
Image of page 6
Image of page 7

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

View Full Document Right Arrow Icon
Image of page 8
Image of page 9

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

View Full Document Right Arrow Icon
Image of page 10
Image of page 11

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

View Full Document Right Arrow Icon
Image of page 12
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: ”Ofe‘é £v Tim's The University of Queensland School of Engineering MECH2700 Engineering Analysis I (2005) Exercise Sheet 7 1; A polynomial passes through the points (3.0,2.0), (4.0,3.1), (5.0,4.4) and (6.0,6.0) in an (11:, y) coordinate system. The polynomial can be written as y=ao+a1x+a2x2+a3x3 Set up the system of linear equations which specifies that the polynomial passes through the points and solve for the coefficients (1.0 through a3. . The matrix for this problem is called a Vandermonde matrix. Larger matrices of this form tend to be ill conditioned. What is the condition number of the 4 x 4 matrix for the set of 4 interplation points? If two more points (7.0, 7.0) and (8.0, 8.0) are added to the set, what is the condition number of the new (6 X 6) matrix. Plot the two sets of points and their interpolating polynomials. 2. Write a Python function that uses the Lagrange basis to interpolate a set of (x, y) coordinates. Your function should accept a list of coordinates and the value of x at which you want the value of y from your interpolant. Be sure to document your function sothat others who use it will know what is expected and what will be returned. 3. Write a Python function to integrate user-defined functions f (:13) over specified ranges of m using the trapezoidal rule with n panels. A suitable prototype might be def trap_integrate(f; a, b, n): # ... # ~—— guts of your function here --- # . . . return integrand where integrand is an approximation to the integral (2 f f (x) dzz: . Try your integration function on the following "sinrc 7', I12] —-—dac, 12:3] smmdx o 56 0 You will need to handle the a: = 0 point carefully for the first function. How large a value for n isrequired to get an accuracy of 4 significant figures? 1 4. Rework the previous exercise with the Newton-Cotes 5-point formula instead of the trapeziodal rule. {:74 2h , 8 7 (6) / f(x)drr = —— (7f0 + 32f1 + 12 f2 + 32 f3 + m) + ——h f (5) mo 45 945 . Write a Python function to integrate f ((1)) over the range [a, b] and have this function adaptively subdivide the original range to control the size of the truncation error. It turns out to be relatively easy to use recursive ftmction calls to achieve this. A rough estimate of the truncation error of the integrand on each segment can be made by comparing the five-point estimate (given above) with the three-point estimate 1 5 (4) 901» f (6) f: f(z)d:v = 5;- (fo + 4f1 + f2) + Since we already have to evaluated the functions at these points for the 5—point estimate, this isn’t much extra work. If the difference is smaller than a certain tolerance, return the 5—point estimate, else split the segment into two equal parts and apply the function to each of the parts, returning the sum of the two. Apply your adaptive integration function to the two integrals Il=f01fii2§da 12:]: 1 dm 1+:z:2 with a tolerance of 1.0 X 10‘6 at each segment. How many function evaluations are required and what are the final errors. (Note that each of the integrals evaluates to 7r/4.) # file: poly_fit.py # MECH27OO Ex sheet 7, Q l. # PJ, 20—may—04 from.poly import horner from.matrix import * print “First set of data points..." xa = [3.0, 4.0, 5.0, 6.0] ya = [2.0, 3.1, 4.4, 6.0] n = len(xa) vdm = new_matrix(n, n) for i in range(n): for j in range(n): vdei][j] = xa[i]**j b = transpose([ya]); # a vector is represented as a list of a single list alpha = gauss~jordan(vdm, b) # get our coefficients back into a single list for later use alist = U for i in range(n): alist.append(alpha[i][0]) print “coefficients=", alist print “condition number=‘, cond(vdm) # Sample data for later plotting. xpl = [1;yp1 = [1 dx = 0.1 . for i in range(81): x = 2.0 + dx * i y = horner(x, alist) xpl.append(x); yp1.append(y) print “Add two data points...“ Xa.append(7.0); ya.append(7.0) xa.append(8.0); Ya.append(8.0) n = len(xa) vdm = new;matrix(n, n) for i in range(n): for j in range(n): vdm[i][j] = xa[i]**j b = transpose([ya]); # a vector is represented as a list of a single list alpha = gauss_jordan(vdm, b) # get our coefficients back into a single list for later use alist = [J for i in range(n): alist.append(alpha[i][0]) print “coefficients=", alist print “condition number=", cond(vdm) # Sample data for later plotting. XpZ = [1; yp2 = [1 dx = 0.1 for i in range(81): x = 2.0 + dx * i y = horner(x, alist) xp2.append(x); yp2.append(y) # _-_._ __ ____ _ _ # Create the plot of both fitted polynomials and the original data import Gnuplot grph = Gnuplot.Gnuplot() dpoints = Gnuplot.Data(xa, ya, with=’points’) d1 = Gnuplot.Data(xpl, ypl, with=’lines’) d2 = Gnuplot.Data(xp2, yp2, with=’lines') grph.title(’ExerciSe Sheet 7, Q1, polynomial fit') grph(’set yrange [—5:20]’) grph.xlabel('x'); grph.ylabel(’y') grph.plot(dpoints, d1, d2) grph.hardcopy(’poly_fit.eps’, terminal=’postscript’, mode='eps', fontsize=20) First set of data points... coefficients: [—1.0999999999999857, 1.1833333333333222, —0.099999999999997202, 0.0166666666666664 41] condition number: 28749.0 Add two data points... coefficients= [—68.600000000014404, 72.808333333348244, —29.537500000005991, 5.870833333334482, — 0.56250000000010525, 0.020833333333337034] condition number= 38986905.6 Exercise Sheet 7, Q1, polynomial fit # POlY-PY def polyn(x, a): “Returns p(x)" p = 0.0 for i in range(len(a)): p += a[i] * x**i return p def horner(x, a): "Returns p(x) evaluated via Horner’s rule." n = len(a) indx = range(n~l) indx.reverse() p = a[n—l] for i in indx: p = a[i] + x * p return p # ____ _ .. _. _ _ _ _ if __name__ == ’__main__’: print "Start poly..." result = polyn(0.l, [1.0, 2.0, 3.0, 4.0, 5.0]) print "result=“, result print "Horner’s rule" result2 = horner(0.l, [1.0, 2.0, 3.0, 4.0, 5.0]) print "result2=", resultZ print "end." # file: lagrange.py "I'll Denmnstration code for Lagrange interpolation. PJ, MECH27oo, May 2004 Transcript... Start Lagrange... x= 2.4 result= 1.28 end "'1" def lbasistj, x, xa): "Returns Lagrange basis fn j n = len(xa) p = 1.0 for i in range(n): if i <> j: p *= return p def lpoly(x, xa, ya): “Returns the value of the Lagrange polynomial for data points xa,ya." sum = 0.0 for j in range(len(xa)): (x — xaEiD / (xa[j] sum_+= ya[j] * lbasis(j, x, return sum if __name~_ == ’__main__’: print ”Start Lagrange..." Xx = [1.0, 2.0, 3.0] yy = [1.0, 1.0, 2.0] x = 2.4 result = lpoly(x,xx,yy) print "x=", x, "result=", print "end" result (for array xa) evaluated at x" — Xa[i]) xa) N M 5m) =1.85193¥o (Tabb. 5‘3 Abramwfii % s’f‘egun ) Q4 MM ham :1 5fm‘dm Wang/med“ Cymbal/Mumgmnmjg, # file: trap_integrate.py "II" a. Quadrature using the trapezoidal rule —~ Ex Shit 7, Q3. PJ, 2 O—May—O 4 Expected output... n 11 errorl 12 error2 10 1.849317 —2.620e—03 1.983524 —l.648e—02 129 get we M451. W“ 50 1.851832 —1.047e—O4 1.999342 ~6.580e—O4 9/ 100 1.851911 —2.613e—05 1.999836 —1.645e-—O4 150 1 851925 —l.158e—O5 1.999927 —7.311e—05 200 1.851931 —6.493e-06 1.999959 —4.112e—05 400 1.851935 —1.584e—06 1.999990 ~1.028e—05 600 1.851936 —6.752e~07 1.999995 ~4.569e—06 800 1.851937 —3.571e—07 1.999997 —2.570e—06 1000 1.851937 —2.098e—07 1.999998 —1.645e—06 def trap_integrate(f, a, b, n): "Returns the integral of f(x) from.x=a to x=b using n panels.“ assert callable(f) assert n >= 1 h = (b — a) / n integrand = 0.0 x = a; fiml = f(x) for i in range(l,n+l): X = a + h * i; fi = f(x) integrand += 0.5 * h * (fiml + fi) fiml = fi # save for next iteration return integrand if __name__ == ’«_main ': import math def f1(x): if x != 0.0: return math.sin(x) / x else: return 1,0 def f2(x): return math.sin(x) exactl = 1.8519370 # from table 5.3 Abramowitz and Stegun exath = 2.0 print “ n 11 errorl I2 errorZ" print " ————————————————————————————————————————————————— " for n in [10, 50, 100, 150, 200, 400, 600, 800, 1000]: I1 = trap_integrate(fl, 0.0, math.pi, n) 12 = trap_integrate(f2, 0.0, math.pi, n) print “%4d %10.6f %10.3e %10.6f %10.3e“ % \ (n, 11, Il-exactl, IZ, IZ—exath) print " — —————+ ~—— — ——————————————————————————— " print “Done." # f ile: ncote35_integrate.py "I!" Quadrature using the Newton~Cotes 5—point rule -— Ex Sheet 7, Q4. PJ, 20—Mayf04 Expected output... def def n 11 errorl 12 error2 1 1.851897 —4.027e—05 1.998571 ~1.429e—03 2 1.851937 —4.820e—07 l.999983 ~l.687e—05 3 1.851937 6.4llen09 1.999999 —1.413e—06 4 1.851937 4.395e—08 2.000000 —2.475e—07 6 WW? ‘Ff/w segmomfis OJZL vuaecL;J lgr-ixé Ctfivedzkfi frac&a\$n e. ncotes5_panel(f, a, b): "Returns the integral of f(x) from.x=a to x=b using N—C 5—point rule.“ dx = b — a f0 = f(a); fl = fta + 0.25*dx) f2 = f(a + 0.5*dx) f3 = f(a + 0.75*dx) f4 = f(b) I = (b — a) / 90.0 * (7.0*f0 + 32.0*f1 + 12.0*f2 + 32.0*f3 + 7.0*f4) return I - , ncote55_integrate(f, a, b, n): "Returns the integral of f(x) from.x=a to x=b using n panels.“ assert callable(f) assert n >= 1 h = (b — a) / n integrand = 0.0 for i in range(n): x1 = a + h * 1 x2 = x1 + h integrand += ncotesSkpane1(f, x1, x2) return integrand if name__ == ’ main ': import math def f1(x): if x != 0.0: return math.sin(x) / x else: return 1.0 and Stegun def f2(x): return math.sin(x) exactl = 1.8519370 # from table 5.3 Abramowitz exath = 2.0 print " n I1 errorl IZ errorZ“ print " ——————————————————————————— — ~ ——‘— ———- ——— for n in [1, 2, 3, 4]: I1 = ncotes5_integrate(fl, 0.0, math.pi, n) I2 = ncote55_integrate(f2, 0.0, math.pi, n) print “%4d %10.6f %10.3e %10.6f %10.3e" % \ (n, 11, Il—exactl, I2, IZ—exath) print “ ——— __ — _ _ _ ____1 -1- _ ~_ print “Done." # file: adapti.py “I!“ Adaptive quadrature using Newton—Cotes 5— and 3—point rules. P. Jacobs School of Engineering, UQ 18—Nov-2003 I'll" def rinteg """Apply Newton—Cotes 5— and 3—point quadrature rules to the segment [a,b]. Input. f : a, b: tol : ( f, a, b, tol ): user—supplied function, f(x) range of integration Returns integral of f(x) from a to b. dx = b f0 = f f1 = f f2 = f f3 = f f4 = f 12 = I4 = if abs " a (a) (a + 0.25 * dx) (a + 0.5 * dx) (a + 0.75 * dx) (b) dX/6.0 * (f0 + 4 * £2 + £4) dx/90.0 * (7*f0 + 32*f1 + 12*f2 + 32*f3 + 7*f4) (I4 ~ I2) > tol: mid=0.5 * (a+b) I else: I = I4 return I from math import pi, sqrt if name == “ main “: countl = 0 def funl(x): g1 obal countl countl += 1 if abs(x) < 1.0: return sqrt(l.0 — x * x) else: count2 return 0.0 = 0 def fun2(x): global count2 countZ += 1 return 1.0 / (1.0 + x * x) print a = O. pi4_l pi4_2 print print print print “Begin adapti...“ 0; b = 1.0; = rinteg(funl, a, b, 1.0e-6) = rinteg(fun2, a, b, 1.0e—6) "Estimates of pi/4: ", pi4_l, p14_2 "errors:“, pi/4 — pi4_1, pi/4 — pi4_2 "number of function calls:“, "Done.“ countl, maximum difference between rules above which the range is = rinteg(f, a, mid, tol) + rinteg(f, mid, b, tol) countZ 'split Begin adapti... Estimates of pi/4: 0.785397761493 0.785398163427 errors: 4.01904837855e—07 -2.95610202983e—ll number of function calls: 165 75 Done. Mme“- (Lott; formula” (‘Abmmwnki 3; Stew ML see) 7‘: - ' . 51¢, ’FCX) dX. =' ~b2: (”Poi-fig. B - 1‘12: kg 442%; ) TMPL'E‘OIHQQ 51215 I h ‘ i S (4) - I ' 10 (7Q dvc = 3 (£0 + 44' + £2) \ C35 1,, { (I?) . 8mm» /3 5 fix) du. = 2;), (150+ 3€,+ $185185) ~ $5 In“; 15(4)<E;) Smoaowa/E” ‘ ) d : &h - 7 __ "l (6) in at 1 ZS— (7%1-32‘E + libel +3>2xc5+ (:4) CST—ESL} 19 . 6;) (such) dx .= $23. (HQ +7543| +SO¥Z+5O€3+751C4+|9155) “Lo 9,68 _ 275 W {3&2}? ) \ZOGJG x5» - - . 3; (70013:, = L (4\£+ alga + 37151 + 2733+ a7 a + 244, g +4HEB> ° . 140 V _ 91,19 “C(23). 1400 (2) 3C7 ., - ‘ v - S «Cm dx = 7h- ‘WSlQCOHCfl + 3577 (£446) +132?>(4’2+{S 35+ 2ng (($104 )1! 1° N 2.80 ‘ \ 818% W ya) (2) 5V84—Oo L3 . . S (la) d1..— 4*: {chx-xfla)+5888(fi+€7)~925(¥2+€6) 1° l4r15‘ - + \otmlc (€34?ch ‘ 4540 {or} ‘ (if-1:, (my?) (Note, *Hnai“ 'HNL Q‘PDNC' 4’0)”me has: 50mm. mga‘jfia ma‘gkts,) ...
View Full Document

{[ snackBarMessage ]}

What students are saying

  • Left Quote Icon

    As a current student on this bumpy collegiate pathway, I stumbled upon Course Hero, where I can find study resources for nearly all my courses, get online help from tutors 24/7, and even share my old projects, papers, and lecture notes with other students.

    Student Picture

    Kiran Temple University Fox School of Business ‘17, Course Hero Intern

  • Left Quote Icon

    I cannot even describe how much Course Hero helped me this summer. It’s truly become something I can always rely on and help me. In the end, I was not only able to survive summer classes, but I was able to thrive thanks to Course Hero.

    Student Picture

    Dana University of Pennsylvania ‘17, Course Hero Intern

  • Left Quote Icon

    The ability to access any university’s resources through Course Hero proved invaluable in my case. I was behind on Tulane coursework and actually used UCLA’s materials to help me move forward and get everything together on time.

    Student Picture

    Jill Tulane University ‘16, Course Hero Intern