Computational Turbulent Reacting Flow

330T Engineering Building
Brigham Young University
Provo, UT 84602
801-422-1772
davidlignell@byu.edu

Numerical Integration 1


  • Integrate functions in 1D
  • General problem. Find $I$:

$$I = \int_{x_1}^{x_2}f(x)dx.$$

  • ODE solvers do this:

$$\frac{dy}{dx} = f(x) \rightarrow y_2 - y_1 = \int_{x_1}^{x_2}f(x)dx.$$

  • Numerical integration methods have some similiarities to ODE solvers.

Outline

  • Direct fit polynomial
  • Trapezoid rule
  • Simpson’s rule

Direct fit polynomial

  • This is useful when we have an unknown function, but whose values we know at specific given points.
  • Fit a polynomial to the function, then integrate the polynomial between the desired bounds.
    • The fitting part is tricky.
    • The integrating part is easy (analytic).
    • Be wary of high order polynomials.
      • If you fit a high order polynomia to a collection of points, the polynomial may vary widely between the given points, and diverge quickly before and after the first and last points.
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
%matplotlib inline

xgiven = np.arange(11)
ygiven = np.array([1,2,3,1,4,6,1,5,10,2,2])

p = np.polyfit(xgiven,ygiven,len(xgiven)-1)
f = interp1d(xgiven, ygiven)

x       = np.linspace(0,10,1000)
ypoly   = np.polyval(p,x)
yinterp = f(x)

fig, (ax1, ax2, ax3) = plt.subplots(3,1, figsize=(5,10))

ax1.fill_between(x, 0, yinterp, facecolor='lightgray')
ax1.plot(x,yinterp, 'r-')
ax1.plot(xgiven,ygiven, 'ko')
ax1.set_ylim([0,13])
ax1.set_title('linear interpolation')
ax1.tick_params(labelbottom=False)    

ax2.fill_between(x, 0, ypoly, facecolor='lightgray')
ax2.plot(x,ypoly, 'b-')
ax2.plot(xgiven,ygiven, 'ko')
ax2.set_ylim([0,13])
ax2.set_title('polynomial fit')
ax2.tick_params(labelbottom=False)    

ax3.plot(x,yinterp, 'r-')
ax3.plot(x,ypoly, 'b-')
ax3.fill_between(x, yinterp, ypoly, facecolor='lightgray')
ax3.plot(xgiven,ygiven, 'ko')
ax3.set_ylim([0,13])
ax3.set_title('difference');

Trapezoid Rule

  • Use a uniform grid of spacing $\Delta x$.
  • Approximate the function as linear on each interval.
x    = np.arange(6)
fpts = np.array([2, 1, 2.5, 3.5, 1, 1.75])
Fi = interp1d(x, fpts)
Fs = interp1d(x, fpts, kind='cubic')
xx = np.linspace(0,5,1000)
fi = Fi(xx)
fs = Fs(xx)

plt.rc('font', size=12)
fig, ax = plt.subplots(1,1)
ax.set_ylim([0,4.1])
ax.plot(xx,fs, 'k-', label='function')
ax.plot(xx,fi, 'r-', label='trapezoid', linewidth=0.75)
ax.fill_between(xx,0,fi, facecolor='lightgray')
ml, sl, bl = ax.stem(x,fpts, bottom=-1)
plt.setp(sl, 'color', 'k')
plt.setp(ml, 'color', 'k')
ax.legend(frameon=False)

for i in x[:-1]:
    ax.text(x[i]+0.4,0.5,r'$I_'+str(x[i])+'$')
for i in x[:]:
    ax.text(x[i]-0.1,fpts[i]+0.3,r'$f_'+str(x[i])+'$')
  • Now we have

$$I = I_0 + I_1 + I_2 + I_3 + I_4.$$

$$I = -\frac{\Delta x}{2}(f_0+f_{n-1}) + \Delta x\sum_{i=0}^{n-1}f_i.$$
  • Note the overlap of the terms.
    • There are two of every function value except the first and the last.
    • So, we double all of them and then subtract the first and last.
  • This method can be used for nonuniform grids (but then we can’t simply factor out the $\Delta x$.
  • The method is 2$^{nd}$ order. That is, the global error is $\propto \Delta x^2$.

Midpoint method

  • Another second order method is the midpoint method.
  • The trapezoid method effectively takes the sum of the areas of intervals, where each interval area is the width times the function evaluated as the average of its values at the endpoints of the intervals.
  • The midpoint method takes the sum of the areas of intervals, where each interval area is the width times the function evaluated at the midpoint of the interval.

Simpson’s Rule

  • Instead of fitting a line between adjacent points, as in the Trapezoid Rule, we fit a parabola through 3 adjacent points.
n    = 7
x    = np.arange(n)
fpts = np.array([2, 1, 2.5, 3.5, 1, 1.75, 3])
dx   = x[1]-x[0]
Fs = interp1d(x, fpts, kind='cubic')
xx = np.linspace(0,n-1,1000)
fs = Fs(xx)

plt.rc('font', size=12)
fig, ax = plt.subplots(1,1)
ax.set_ylim([0,5.0])
ax.plot(xx,fs, 'k-')

nI = int((n-1)/2)
a = np.zeros(nI); b=a.copy(); c=a.copy()
for i in range(0,nI):
    a[i] = 1/2/dx**2*(fpts[0+2*i] -2*fpts[1+2*i] +     fpts[2+2*i])
    b[i] = 1/dx*(-1.5*fpts[0+2*i] +2*fpts[1+2*i] - 0.5*fpts[2+2*i])
    c[i] = fpts[0+i*2]
    xi   = np.linspace(i*2,(i+1)*2,1000)
    fi   = a[i]*(xi-i*2)**2 + b[i]*(xi-i*2) + c[i]
    ax.plot(xi,fi, 'r-', linewidth=0.75)
    ax.fill_between(xi,0,fi, facecolor='lightgray')
    
ml, sl, bl = ax.stem(x[[0,2,4,6]],fpts[[0,2,4,6]], bottom=-1)
ax.plot(x,fpts,'ko')
plt.setp(sl, 'color', 'k')
plt.setp(ml, 'color', 'k')
ax.legend(['function', 'simpson'], frameon=False)

for i in x[[0,2,4]]:
    ax.text(x[i]+0.9,0.5,r'$I_'+str(x[i])+'$')
for i in x[:]:
    ax.text(x[i]-0.1,fpts[i]+0.3,r'$f_'+str(x[i])+'$')
  • Consider each segment individually.
  • Then $x_0=0$, $x_1=\Delta x$, and $x_2 = 2\Delta x$.
  • Parabola:

$$f(x) = ax^2 + bx + c.$$

  • Set $f(0)=f_0$, $f(\Delta x)=f_1$, $f(2\Delta x)=f_2$. This gives three equations for each of our unknowns $a$, $b$, and $c$:

  • Solving these yeilds:

  • Now, the integral over the first interval is:

  • Finally, insert $a$, $b$, and $c$ into $I_0$, and simplify:

$$I_0 = \frac{\Delta x}{3}(f_0 + 4f_1 + f_2).$$

  • Now, we can repeat this process for the other segments and then add all intervals:

$$I = I_0 + I_1 + I_2.$$

$$I = \frac{\Delta x}{3}\left[\underbrace{f_0 + 4f_1 + f_2}{\text{segment 1}} + \underbrace{f_2 + 4f_3 + f_4}{\text{segment 2}} + \underbrace{f_4 + 4f_5 + f_6}_{\text{segment 3}}\right].$$

  • Note, between segments we can combine the recurring terms: $f_2+f_2\rightarrow 2f_2$, and $f_4 + f_4\rightarrow 2f_4$.
  • The coefficients of $f$ are either 1, 4, or 2.

$$I = \frac{\Delta x}{3}\left[f_0 + 4f_1 + 2f_2 + 4f_3 + 2f_4 + 4f_5 + f_6\right].$$

  • Group terms:
$$I = \frac{4\Delta x}{3}\sum_{i\,\text{odd}}f_i + \frac{2\Delta x}{3}\sum_{i\,\text{even}}f_i - \frac{\Delta x}{3}(f_0 + f_{n-1}).$$

  • Here, we have done $\pm\frac{\Delta x}{3}f_0$ and $\pm\frac{\Delta x}{3}f_{n-1}$ to get 2 as the coefficient of $f_0$ and $f_{n-1}$ so that the even sum can go over all even points; this then results in the final correction term.

  • Simpson’s rule requires an odd number of points: $n$ is odd.

Question: what if we have an even number of points?

  • Could do Simpson’s rule for $n-1$ points and Trapazoid rule for the last interval.

    • But that is lower order.

  • Could do Simpson’s rule for $n-1$ points, then integrate a parabola through the last three points over the last sub interval.

  • Could use Simpson’s 3/8 rule.

    • That’s a 4 point method with the same order. So apply that to the first 4 points, then the usual Simpson’s rule to the rest of the intervals.

    $$\int_{x_0}^{x_3}f(x)dx = \frac{3}{8}\Delta x(f_0 + 3f_1 + 3f_2 + f_3).$$

Question: what if the endpoints are hard?

  • For example, if the function goes to $-\infty$ as $x\rightarrow x_0$.
  • Can extrapolate the function.
  • Assume a linear relation in the region $[x_0,,x_1]$ using the function between $[x_1,,x_2]$.

Hence,

$$I_0 = \int_{x_0}^{x_1}f(x)dx = \frac{a\Delta x^2}{2} + b\Delta x,$$

Insert $a$ and $b$ to get:

$$I_0 = \Delta x\left(\frac{3}{2}f_1 - \frac{1}{2}f_2\right).$$

  • For a parabola, computed as for Simpson’s, we would get

$$I_0 = \Delta x\left(\frac{23}{12}f_1 - \frac{16}{12}f_2 + \frac{5}{12}f_3\right).$$