# Simpson's integration method (based on the application "Numerical integrals with various rules" by Marshall Hampton and Nick Alexander)
# Lauri Ruotsalainen, 2010
def parabola(a, b, c):
A, B, C, x = var("A, B, C, x")
K = solve([A*a[0]^2+B*a[0]+C==a[1], A*b[0]^2+B*b[0]+C==b[1], A*c[0]^2+B*c[0]+C==c[1]], [A, B, C], solution_dict=True)[0]
f = K[A]*x^2+K[B]*x+K[C]
return f
@interact
def simpson(
f = input_box(default = x*sin(x)+x+1),
n=slider(2,100,2,6),
interval=input_box(default=(0,10), label="Integrointiväli")
):
f(x) = f
xs = []; ys = []
dx = (interval[1]-interval[0])/n
for i in range(n+1):
xs.append(interval[0] + i*dx)
ys.append(f(xs[-1]))
parabolas = Graphics()
lines = Graphics()
for i in range(0, n-1, 2):
p(x) = parabola((xs[i],ys[i]),(xs[i+1],ys[i+1]),(xs[i+2],ys[i+2]))
parabolas += plot(p(x), x, xmin=xs[i], xmax=xs[i+2], color="red")
lines += line([(xs[i],ys[i]), (xs[i],0), (xs[i+2],0)],color="red")
lines += line([(xs[i+1],ys[i+1]), (xs[i+1],0)], linestyle="-.", color="red")
lines += line([(xs[-1],ys[-1]), (xs[-1],0)], color="red")
show(plot(f(x),x,interval[0],interval[1]) + parabolas + lines, xmin = interval[0], xmax = interval[1])
numeric_value = integral_numerical(f,interval[0],interval[1])[0]
approx = dx/3 *(ys[0] + sum([4*ys[i] for i in range(1,n,2)]) + sum([2*ys[i] for i in range(2,n,2)]) + ys[n])
sum_formula_html = "\\frac{d}{3} \cdot \\left[ f(x_0) + %s + f(x_{%s})\\right]" % (
join([ "%s \cdot f(x_{%s})" %(i%2*(-2)+4, i+1) for i in range(0,n-1)], " + "),
n
)
sum_placement_html = "\\frac{%s}{3} \cdot \\left[ f(%s) + %s + f(%s)\\right]" % (
dx,
N(xs[0],digits=5),
join([ "%s \cdot f(%s)" %(i%2*(-2)+4, N(xk, digits=5)) for i, xk in enumerate(xs[1:-1])], " + "),
N(xs[n],digits=5)
)
sum_values_html = "\\frac{%s}{3} \cdot \\left[ %s %s %s\\right]" %(
dx,
"%s + "%N(ys[0],digits=5),
join([ "%s \cdot %s" %(i%2*(-2)+4, N(yk, digits=5)) for i, yk in enumerate(ys[1:-1])], " + "),
" + %s"%N(ys[n],digits=5)
)
html(r'''
<div class="math">
\begin{align*}
\int_{%s}^{%s} {f(x) \, dx} & = %s \\\
\\\
\int_{%s}^{%s} {f(x) \, dx}
& \approx %s \\\
& = %s \\\
& = %s \\\
& = %s
\end{align*}
</div>
''' % (
interval[0],interval[1],
N(numeric_value,digits=7),
interval[0], interval[1],
sum_formula_html, sum_placement_html, sum_values_html,
N(approx,digits=7)
))
