Simpson Iteration -- Sage

We are integrating a function f over an interval [a,b],
we need the maximum m of the fourth derivative of f in that interval \max_{\xi \in [a,b]} {\left| f^{(4)}( \xi )\right|},
eps is the precision

def h4(d): (a,b,M,eps) = d return 180*eps / (M*(b-a)) def h(d): return sqrt(sqrt(h4(d))) def N(d): (a,b,M,eps) = d return (b-a) / h(d) def eps(h4, a, b, M): return h4 * (b-a) * M / 180

Tabulation of some calculations

def row(a,b,M,e): d = (a,b,M,e) if not eps(h4(d),a,b,M) == e: print "this shouldn't happen" return (e, float(h4(d)), float(h(d)), float(N(d)), ceil(N(d))) def table(a,b,M,es): rows = [("eps", "h^4","h","N","ceil(N)")] rows += map(lambda e: row(a,b,M,e), es) return rows

HTML <table> rendering of tabular data

# surround list elements with tags and return single string def tagmap(tag, xs): return ''.join(map(lambda s: "<%s>%s</%s>"%(tag,s,tag), xs)) # build table from rows, with first row as table headers def htmltable(rows): t = [""" <style>table.my, table.my td, table.my th { border: 1px solid gray; border-collapse: collapse; padding: 6px; }</style> <table class='my'>"""] t.append("<tr>%s</tr>" % tagmap("th", rows[0])) for r in rows[1:]: t.append("<tr>%s</tr>" % tagmap("td", r)) t.append("</table>") return html('\n'.join(t))

epsilons = map(lambda i: 1/(10**i), range(12)) htmltable(table(0,1,100,epsilons))




eps h^4 h N ceil(N)
1 1.8 1.15829218529 0.86334002137 1
1/10 0.18 0.651355562433 1.53525978387 2
1/100 0.018 0.366284150148 2.73012086271 3
1/1000 0.0018 0.205976714391 4.85491771707 5
1/10000 0.00018 0.115829218529 8.6334002137 9
1/100000 1.8e-05 0.0651355562433 15.3525978387 16
1/1000000 1.8e-06 0.0366284150148 27.3012086271 28
1/10000000 1.8e-07 0.0205976714391 48.5491771707 49
1/100000000 1.8e-08 0.0115829218529 86.334002137 87
1/1000000000 1.8e-09 0.00651355562433 153.525978387 154
1/10000000000 1.8e-10 0.00366284150148 273.012086271 274
1/100000000000 1.8e-11 0.00205976714391 485.491771707 486




eps h^4 h N ceil(N)
1 1.8 1.15829218529 0.86334002137 1
1/10 0.18 0.651355562433 1.53525978387 2
1/100 0.018 0.366284150148 2.73012086271 3
1/1000 0.0018 0.205976714391 4.85491771707 5
1/10000 0.00018 0.115829218529 8.6334002137 9
1/100000 1.8e-05 0.0651355562433 15.3525978387 16
1/1000000 1.8e-06 0.0366284150148 27.3012086271 28
1/10000000 1.8e-07 0.0205976714391 48.5491771707 49
1/100000000 1.8e-08 0.0115829218529 86.334002137 87
1/1000000000 1.8e-09 0.00651355562433 153.525978387 154
1/10000000000 1.8e-10 0.00366284150148 273.012086271 274
1/100000000000 1.8e-11 0.00205976714391 485.491771707 486

htmltable(table(0,1,10,epsilons))




eps h^4 h N ceil(N)
1 18.0 2.05976714391 0.485491771707 1
1/10 1.8 1.15829218529 0.86334002137 1
1/100 0.18 0.651355562433 1.53525978387 2
1/1000 0.018 0.366284150148 2.73012086271 3
1/10000 0.0018 0.205976714391 4.85491771707 5
1/100000 0.00018 0.115829218529 8.6334002137 9
1/1000000 1.8e-05 0.0651355562433 15.3525978387 16
1/10000000 1.8e-06 0.0366284150148 27.3012086271 28
1/100000000 1.8e-07 0.0205976714391 48.5491771707 49
1/1000000000 1.8e-08 0.0115829218529 86.334002137 87
1/10000000000 1.8e-09 0.00651355562433 153.525978387 154
1/100000000000 1.8e-10 0.00366284150148 273.012086271 274




eps h^4 h N ceil(N)
1 18.0 2.05976714391 0.485491771707 1
1/10 1.8 1.15829218529 0.86334002137 1
1/100 0.18 0.651355562433 1.53525978387 2
1/1000 0.018 0.366284150148 2.73012086271 3
1/10000 0.0018 0.205976714391 4.85491771707 5
1/100000 0.00018 0.115829218529 8.6334002137 9
1/1000000 1.8e-05 0.0651355562433 15.3525978387 16
1/10000000 1.8e-06 0.0366284150148 27.3012086271 28
1/100000000 1.8e-07 0.0205976714391 48.5491771707 49
1/1000000000 1.8e-08 0.0115829218529 86.334002137 87
1/10000000000 1.8e-09 0.00651355562433 153.525978387 154
1/100000000000 1.8e-10 0.00366284150148 273.012086271 274

htmltable(table(0,1,1,epsilons))




eps h^4 h N ceil(N)
1 180.0 3.66284150148 0.273012086271 1
1/10 18.0 2.05976714391 0.485491771707 1
1/100 1.8 1.15829218529 0.86334002137 1
1/1000 0.18 0.651355562433 1.53525978387 2
1/10000 0.018 0.366284150148 2.73012086271 3
1/100000 0.0018 0.205976714391 4.85491771707 5
1/1000000 0.00018 0.115829218529 8.6334002137 9
1/10000000 1.8e-05 0.0651355562433 15.3525978387 16
1/100000000 1.8e-06 0.0366284150148 27.3012086271 28
1/1000000000 1.8e-07 0.0205976714391 48.5491771707 49
1/10000000000 1.8e-08 0.0115829218529 86.334002137 87
1/100000000000 1.8e-09 0.00651355562433 153.525978387 154




eps h^4 h N ceil(N)
1 180.0 3.66284150148 0.273012086271 1
1/10 18.0 2.05976714391 0.485491771707 1
1/100 1.8 1.15829218529 0.86334002137 1
1/1000 0.18 0.651355562433 1.53525978387 2
1/10000 0.018 0.366284150148 2.73012086271 3
1/100000 0.0018 0.205976714391 4.85491771707 5
1/1000000 0.00018 0.115829218529 8.6334002137 9
1/10000000 1.8e-05 0.0651355562433 15.3525978387 16
1/100000000 1.8e-06 0.0366284150148 27.3012086271 28
1/1000000000 1.8e-07 0.0205976714391 48.5491771707 49
1/10000000000 1.8e-08 0.0115829218529 86.334002137 87
1/100000000000 1.8e-09 0.00651355562433 153.525978387 154

Simpson Iteration

124 days ago by wires