Sovellus: Simpsonin menetelmä

# Sovellus: Simpsonin menetelmä # Lauri Ruotsalainen, 2011 (perustuu sovellukseen "Numerical integrals with various rules" by Marshall Hampton and Nick Alexander) # Sellaisen paraabelin funktio, joka lukee pisteiden a, b ja c kautta. def paraabeli(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), vali = input_box(default=(0,10), label="Integrointiväli")): f(x) = f xs = []; ys = [] dx = (vali[1]-vali[0])/n # Piirretään paraabelit. for i in range(n+1): xs.append(vali[0] + i*dx) ys.append(f(xs[-1])) paraabelit = Graphics() viivat = Graphics() for i in range(0, n-1, 2): p(x) = paraabeli((xs[i], ys[i]), (xs[i+1], ys[i+1]), (xs[i+2], ys[i+2])) paraabelit += plot(p(x), x, xmin=xs[i], xmax=xs[i+2], color="red") viivat += line([(xs[i], ys[i]), (xs[i], 0), (xs[i+2], 0)], color="red") viivat += line([(xs[i+1], ys[i+1]), (xs[i+1], 0)], linestyle="-.", color="red") viivat += line([(xs[-1],ys[-1]), (xs[-1],0)], color="red") # Esitetään funktion käyrä ja paraabelit kuvaajassa. show(plot(f(x), x, vali[0], vali[1]) + paraabelit + viivat, xmin = vali[0], xmax = vali[1]) # Lasketaan määrätyn integraalin approksimaatio Simpsonin menetelmän ja Sagen tarkan numeerisen metodin avulla. numeerinen_arvo = integral_numerical(f,vali[0],vali[1])[0] approksimaatio = 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]) # Esitetään approksimaation laskemisen välivaiheet. summa_kaava_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) summa_sijoitus_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)) summa_arvot_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> ''' % (vali[0], vali[1], N(numeerinen_arvo, digits=7), vali[0], vali[1], summa_kaava_html, summa_sijoitus_html, summa_arvot_html, N(approksimaatio, digits=7)))