Simpson's Rule

136 days ago by david.guichard

We want to compute the area under the parabola through three points, (x_0,y_0), (x_1,y_1), (x_2,y_2), where x_0=x_1-h and x_2=x_1+h. First, we get Sage to solve for the coefficients a, b, and c of the parabola q(x)=ax^2+bx+c.
var("a b c y0 y1 y2 x1 h"); q(x)=a*x^2+b*x+c;soln=solve([y0==q(x1-h),y1==q(x1),y2==q(x1+h)],a,b,c,solution_dict=true);soln 
        
\newcommand{\Bold}[1]{\mathbf{#1}}\left(a, b, c, y_{0}, y_{1}, y_{2}, x_{1}, h\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{c : -\frac{{\left(2 \, y_{1} - y_{2}\right)} x_{1}^{2} - 2 \, h^{2} y_{1} + h x_{1} y_{2} - {\left(h x_{1} + x_{1}^{2}\right)} y_{0}}{2 \, h^{2}}, b : \frac{2 \, {\left(2 \, y_{1} - y_{2}\right)} x_{1} - {\left(h + 2 \, x_{1}\right)} y_{0} + h y_{2}}{2 \, h^{2}}, a : \frac{y_{0} - 2 \, y_{1} + y_{2}}{2 \, h^{2}}\right\}\right]
\newcommand{\Bold}[1]{\mathbf{#1}}\left(a, b, c, y_{0}, y_{1}, y_{2}, x_{1}, h\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{c : -\frac{{\left(2 \, y_{1} - y_{2}\right)} x_{1}^{2} - 2 \, h^{2} y_{1} + h x_{1} y_{2} - {\left(h x_{1} + x_{1}^{2}\right)} y_{0}}{2 \, h^{2}}, b : \frac{2 \, {\left(2 \, y_{1} - y_{2}\right)} x_{1} - {\left(h + 2 \, x_{1}\right)} y_{0} + h y_{2}}{2 \, h^{2}}, a : \frac{y_{0} - 2 \, y_{1} + y_{2}}{2 \, h^{2}}\right\}\right]

The coefficients are a bit messy. Now we substitute the values into ax^2+bx+c to get the actual parabola.
p=q.substitute(soln[0]);p 
        
\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ \frac{{\left(y_{0} - 2 \, y_{1} + y_{2}\right)} x^{2}}{2 \, h^{2}} + \frac{{\left(2 \, {\left(2 \, y_{1} - y_{2}\right)} x_{1} - {\left(h + 2 \, x_{1}\right)} y_{0} + h y_{2}\right)} x}{2 \, h^{2}} - \frac{{\left(2 \, y_{1} - y_{2}\right)} x_{1}^{2} - 2 \, h^{2} y_{1} + h x_{1} y_{2} - {\left(h x_{1} + x_{1}^{2}\right)} y_{0}}{2 \, h^{2}}
\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ \frac{{\left(y_{0} - 2 \, y_{1} + y_{2}\right)} x^{2}}{2 \, h^{2}} + \frac{{\left(2 \, {\left(2 \, y_{1} - y_{2}\right)} x_{1} - {\left(h + 2 \, x_{1}\right)} y_{0} + h y_{2}\right)} x}{2 \, h^{2}} - \frac{{\left(2 \, y_{1} - y_{2}\right)} x_{1}^{2} - 2 \, h^{2} y_{1} + h x_{1} y_{2} - {\left(h x_{1} + x_{1}^{2}\right)} y_{0}}{2 \, h^{2}}

Now we integrate from x_0=x_1-h to x_2=x_1+h.
integral(p,x,x1-h,x1+h) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ -\frac{{\left(h^{3} - 3 \, h x_{1}^{2} + 2 \, x_{1}^{3}\right)} y_{2} - 4 \, {\left(2 \, h^{3} - 3 \, h^{2} x_{1} + x_{1}^{3}\right)} y_{1} - {\left(5 \, h^{3} - 3 \, h x_{1}^{2} - 2 \, x_{1}^{3}\right)} y_{0}}{12 \, h^{2}} - \frac{{\left(h^{3} - 3 \, h x_{1}^{2} - 2 \, x_{1}^{3}\right)} y_{0} - 4 \, {\left(2 \, h^{3} + 3 \, h^{2} x_{1} - x_{1}^{3}\right)} y_{1} - {\left(5 \, h^{3} - 3 \, h x_{1}^{2} + 2 \, x_{1}^{3}\right)} y_{2}}{12 \, h^{2}}
\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ -\frac{{\left(h^{3} - 3 \, h x_{1}^{2} + 2 \, x_{1}^{3}\right)} y_{2} - 4 \, {\left(2 \, h^{3} - 3 \, h^{2} x_{1} + x_{1}^{3}\right)} y_{1} - {\left(5 \, h^{3} - 3 \, h x_{1}^{2} - 2 \, x_{1}^{3}\right)} y_{0}}{12 \, h^{2}} - \frac{{\left(h^{3} - 3 \, h x_{1}^{2} - 2 \, x_{1}^{3}\right)} y_{0} - 4 \, {\left(2 \, h^{3} + 3 \, h^{2} x_{1} - x_{1}^{3}\right)} y_{1} - {\left(5 \, h^{3} - 3 \, h x_{1}^{2} + 2 \, x_{1}^{3}\right)} y_{2}}{12 \, h^{2}}

Still messy, but let's tell Sage to multiply everything out and simplify.
expand(_) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ \frac{1}{3} \, h y_{0} + \frac{4}{3} \, h y_{1} + \frac{1}{3} \, h y_{2}
\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ \frac{1}{3} \, h y_{0} + \frac{4}{3} \, h y_{1} + \frac{1}{3} \, h y_{2}

Voila! A simple formula for the area, namely (h/3)(y_0+4y_1+y_2).