# Riemann Sum Optimization

At the moment, I am taking a calculus course at my high school (we are starting to learn about integrals) and thought of an interesting problem while learning about Riemann Sums. Once I develop a method on how to solve this problem for any general curve, I hope to write a script in python.

Problem

Consider the following parameters:

• $y = f(x) = (x - 1.5)^{1/3} + 2$
• Domain: $[-2, 6]$
• $n = 4$
• The width of each subdivision, $\Delta w_k$, does not have to be equal.

Using left-hand endpoints, right-hand endpoints, midpoints, or a combination of all three, how could you orient the rectangles to cover the greatest area underneath the curve? Ideally, the intended python script I want to write would involve an algorithm that can be applied to any particular curve, any specific domain interval, and any number of subdivisions.

How would you even approach this problem? What steps would you outline?

Graph of $y = f(x) = (x - 1.5)^{1/3} + 2$

Same Graph with Varying Widths

This is an example that involves the midpoint approximation with varying widths ($\Delta w_k$).

-
I'm not sure I understood what you mean, but script-wise, you can divide your domain to $N$ parts (ie make a vector from beginning to end with some $dt$ resolution) and then take the value of the function at each point. Multiplied by $dt$ you'll get the value of the rectangle whose width is $dt$ and height is the value of the function. As $dt$ gets close to 0, the sum of all rectangles approximates the integral. –  yoki Jan 17 '13 at 22:34
Unfortunately I require at least a 10 reputation to provide images... When I do get a 10 rep, then I will post some images to show you what I mean. –  Biff Jan 17 '13 at 22:36
the $dt$ that you mentioned would be the equivalent of $\Delta w_k$ in this instance. –  Biff Jan 17 '13 at 22:42
There we go, the graphs have been given. –  Biff Jan 18 '13 at 2:26
I understand the graphs, I didn't understand what you are trying to accomplish. I believe I specified an initial algorithm for the calculation in my first message... –  yoki Jan 18 '13 at 9:11

OK, we know that the Riemann sum is an approximation to the integral $I = \int_a^b f(x) \, dx$. In this case, the value of the integral is $I = 16 + (3/4) ((9/2)^{4/3} - (7/2)^{4/3}) \approx 17.5865$. The Riemann sum using the midpoint rule as you specify is

$$S(x_1,x_2,x_3) = \frac{1}{4} \left [ f \left (\frac{a+x_1}{2} \right )(x_1-a) + f \left (\frac{x_1 + x_2}{2} \right )(x_2-x_1) + f \left (\frac{x_2+x_3}{2} \right )(x_3-x_2) + f \left (\frac{x_3+b}{2} \right )(b-x_3) \right ]$$

We want to minimize the difference between $S$ and $I$. Actually, that's not very accurate: we want to minimize the absolute value of the difference between S and I. And if we have access to a multivariable minimization routine, such as simplex, we could do that. However, if we want to analyze things more analytically for your Python code, we would then choose to minimize the square of the difference between $S$ and I:

$$\Delta(x_1,x_2,x_3) = (I - S(x_1,x_2,x_3))^2$$

Now, if this is your first course in Calculus, then what we do here is going to be new; you typically don't lean these concepts until Calculus III. What we do next is to find the derivatives of $\Delta$ with respect to each of the $x_i$, holding the other $x_j$ constant; this is called a partial derivative with respect to $x_i$. For example:

$$\frac{\partial \Delta}{\partial x_1} = -2 (I - S(x_1,x_2,x_3)) \frac{\partial S}{\partial x_1}$$

$$\frac{\partial S}{\partial x_1} = \frac{1}{4} \left [\frac{1}{2} f' \left (\frac{a+x_1}{2} \right )(x_1-a) + f \left (\frac{a+x_1}{2} \right ) + \frac{1}{2} f' \left (\frac{x_1 + x_2}{2} \right )(x_2-x_1) - f \left (\frac{x_1 + x_2}{2} \right ) \right ]$$

Do the same for the other points $x_2$ and $x_3$. Then set each of the $\frac{\partial \Delta}{\partial x_i} = 0$, which is equivalent to setting $\frac{\partial S}{\partial x_i} = 0$. You'll have 3 equations and 3 unknowns in a nonlinear system of equations. You can try to solve for the $x_i$ analytically, or you can use a numerical routine to find them.

I hope that helps. Good luck!

-