2

I'm currently trying to make a program which will plot a function using matplotlib, graph it, shade the area under the curve between two variables, and use Simpson's 3/8th's rule to calculate the shaded area. However, when trying to print the variable I've assigned to the final value of the integral, it prints a list.

To begin, here's the base of my code:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon

This definition defines the function I will be working with here, a simple polynomial.

def func(x):
    return (x - 3) * (x - 5) * (x - 7) + 85

Here is the function which calculates the area under the curve

def simpson(function, a, b, n):
    """Approximates the definite integral of f from a to b by the
    composite Simpson's rule, using n subintervals (with n even)"""

if n % 2:
    raise ValueError("n must be even (received n=%d)" % n)

h = (b - a) / n #The first section of Simpson's 3/8ths rule
s = function(a) + function(b) #The addition of functions over an interval

for i in range(1, n, 2):
    s += 4 * function(a + i * h)
for i in range(2, n-1, 2):
    s += 2 * function(a + i * h)

return s * h / 3

Now the simpson's rule definition is over, and I define a few variables for simplicity.

a, b = 2, 9  # integral limits
x = np.linspace(0, 10) #Generates 100 points evenly spaced between 0 and 10
y = func(x) #Just defines y to be f(x) so its ez later on

fig, ax = plt.subplots()
plt.plot(x, y, 'r', linewidth=2)
plt.ylim(ymin=0)

final_integral = simpson(lambda x: y, a, b, 100000)

At this point something must have broken down, but I'll include the rest of the code in case you can spot the issue further on.

# Make the shaded region
ix = np.linspace(a, b)
iy = func(ix)
verts = [(a, 0)] + list(zip(ix, iy)) + [(b, 0)]
poly = Polygon(verts, facecolor='0.9', edgecolor='0.5')
ax.add_patch(poly)

plt.text(0.5 * (a + b), 30, r"$\int_a^b f(x)\mathrm{d}x$",
     horizontalalignment='center', fontsize=20)
ax.text(0.25, 135, r"Using Simpson's 3/8ths rule, the area under the curve is: ", fontsize=20)

Here is where the integral value should be printed:

ax.text(0.25, 114, final_integral , fontsize=20)

Here is the rest of the code necessary to plot the graph:

plt.figtext(0.9, 0.05, '$x$')
plt.figtext(0.1, 0.9, '$y$')

ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')

ax.set_xticks((a, b))
ax.set_xticklabels(('$a$', '$b$'))
ax.set_yticks([])

plt.show()

When running this program, you get this graph, and a series of numbers has been printed where the area under the curve should be

Any help here is appreciated. I'm totally stuck. Also, sorry if this is a tad long, it's my first question on the forum.

hakplay
  • 41
  • 1
  • 5
  • Your question is long indeed, but it includes all the necessary parts and a clear problem description. Better a long but complete question than a short incomplete one. So thumps up for this nice question. – ImportanceOfBeingErnest Apr 03 '17 at 21:54

2 Answers2

2

Have you tried feeding your simpson() function the func() directly, as opposed to using the lambda setup?

I think this could work:

final_integral = simpson(func, a, b, 100000)

You might also try:

final_integral = simpson(lambda x: func(x), a, b, 100000)    

What is happening is that y is an array with values func(x), and when you use the expression lambda x: y you are actually creating a constant function of the form f(x) = y = const. Your final_integral is then a list of integrals, where each integrand was the constant function with a particular value from the y array.

Note that you might want to format this number when you print it on the graph, in case it has a lot of trailing decimal points. How you do this depends on whether you are using Python 2 or 3.

ImportanceOfBeingErnest
  • 321,279
  • 53
  • 665
  • 712
1

enter image description hereYou assigned x as linspace which is an array so y is also an array of values of a function of x. You can replace this line of code:

#old:
final_integral = simpson(lambda x:y, a, b, 100000)

#new:
final_integral = simpson(lambda t:func(t), a, b, 100000)

Changing the variable from x to t will give you the value for the area under that the curve. Hope this helps.

jose_bacoy
  • 12,227
  • 1
  • 20
  • 38