

In [1]:

    
from sympy import init_session
init_session()









    



IPython console for SymPy 0.7.7.dev (Python 2.7.8-64-bit) (ground types: gmpy)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/dev

Quadratic

We want to construct a quadratic polynomal through that points $x_{i-1}$, $x_i$, $x_{i+1}$ that gives the correct averages, $f_{i-1}$, $f_i$, and $f_{i-1}$ when integrated over the volume, e.g.

$$\frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} f(x) dx = f_i$$

There are 3 unknowns in the quadratic and three constraints, so this is a linear system we can solve.

Define the quadratic polynomial



In [2]:

    
a, b, c = symbols("a b c")
x0 = symbols("x0")
f = a*(x-x0)**2 + b*(x-x0) + c
f









    Out[2]:





$$a \left(x - x_{0}\right)^{2} + b \left(x - x_{0}\right) + c$$



In [3]:

    
dx = symbols("\Delta")

constraints

Define the 3 constraint equations---here we set them up construct $A$, $B$, and $C$ as the integrals over the 3 control volumes



In [22]:

    
fm, f0, fp = symbols("f_{i-1} f_i f_{i+1}")
#xm32, xm12, xp12, xp32 = symbols("x_{i-3/2} x_{i-1/2} x_{i+1/2} x_{i+3/2}")
xm32 = x0 - Rational(3,2)*dx
xm12 = x0 - Rational(1,2)*dx
xp12 = x0 + Rational(1,2)*dx
xp32 = x0 + Rational(3,2)*dx

interfaces



In [23]:

    
xm32, xm12, xp12, xp32









    Out[23]:





$$\left ( - \frac{3 \Delta}{2} + x_{0}, \quad - \frac{\Delta}{2} + x_{0}, \quad \frac{\Delta}{2} + x_{0}, \quad \frac{3 \Delta}{2} + x_{0}\right )$$



In [5]:

    
A = simplify(integrate(f/dx, (x, xm32, xm12)))
B = simplify(integrate(f/dx, (x, xm12, xp12)))
C = simplify(integrate(f/dx, (x, xp12, xp32)))

The analytic forms of the integrals



In [6]:

    
A, B, C









    Out[6]:





$$\left ( \frac{13 a}{12} \Delta^{2} - \Delta b + c, \quad \frac{\Delta^{2} a}{12} + c, \quad \frac{13 a}{12} \Delta^{2} + \Delta b + c\right )$$

Our linear system is now:

$$A = f_{i-1}$$$$B = f_i$$$$C = f_{i+1}$$

Now find the coefficients of the polynomial



In [7]:

    
coeffs = solve([A-fm, B-f0, C-fp], [a,b,c])
coeffs









    Out[7]:





$$\left \{ a : \frac{1}{2 \Delta^{2}} \left(- 2 f_{i} + f_{{i+1}} + f_{{i-1}}\right), \quad b : \frac{f_{{i+1}} - f_{{i-1}}}{2 \Delta}, \quad c : \frac{13 f_{i}}{12} - \frac{f_{{i+1}}}{24} - \frac{f_{{i-1}}}{24}\right \}$$

And in pretty form, here's the polynomial



In [8]:

    
f.subs(a,coeffs[a]).subs(b,coeffs[b]).subs(c,coeffs[c])









    Out[8]:





$$\frac{13 f_{i}}{12} - \frac{f_{{i+1}}}{24} - \frac{f_{{i-1}}}{24} + \frac{1}{2 \Delta} \left(f_{{i+1}} - f_{{i-1}}\right) \left(x - x_{0}\right) + \frac{\left(x - x_{0}\right)^{2}}{2 \Delta^{2}} \left(- 2 f_{i} + f_{{i+1}} + f_{{i-1}}\right)$$

Cubic

We want to construct a cubic polynomal through that points $x_{i-2}$, $x_{i-1}$, $x_i$, $x_{i+1}$ that gives the correct averages, $f_{i-2}$, $f_{i-1}$, $f_i$, and $f_{i-1}$ when integrated over the volume of each zone



In [10]:

    
a, b, c, d = symbols("a b c d")
f = a*(x-x0)**3 + b*(x-x0)**2 + c*(x-x0) + d
f









    Out[10]:





$$a \left(x - x_{0}\right)^{3} + b \left(x - x_{0}\right)^{2} + c \left(x - x_{0}\right) + d$$

Now perform the integals of $f(x)$ over each zone



In [20]:

    
fm2, fm, f0, fp = symbols("f_{i-2} f_{i-1} f_i f_{i+1}")
xm52 = x0 - Rational(5,2)*dx
xm32 = x0 - Rational(3,2)*dx
xm12 = x0 - Rational(1,2)*dx
xp12 = x0 + Rational(1,2)*dx
xp32 = x0 + Rational(3,2)*dx

interfaces



In [21]:

    
xm52, xm32, xm12, xp12, xp32









    Out[21]:





$$\left ( - \frac{5 \Delta}{2} + x_{0}, \quad - \frac{3 \Delta}{2} + x_{0}, \quad - \frac{\Delta}{2} + x_{0}, \quad \frac{\Delta}{2} + x_{0}, \quad \frac{3 \Delta}{2} + x_{0}\right )$$



In [12]:

    
A = simplify(integrate(f/dx, (x, xm52, xm32)))
B = simplify(integrate(f/dx, (x, xm32, xm12)))
C = simplify(integrate(f/dx, (x, xm12, xp12)))
D = simplify(integrate(f/dx, (x, xp12, xp32)))



In [13]:

    
A, B, C, D









    Out[13]:





$$\left ( - \frac{17 a}{2} \Delta^{3} + \frac{49 b}{12} \Delta^{2} - 2 \Delta c + d, \quad - \frac{5 a}{4} \Delta^{3} + \frac{13 b}{12} \Delta^{2} - \Delta c + d, \quad \frac{\Delta^{2} b}{12} + d, \quad \frac{5 a}{4} \Delta^{3} + \frac{13 b}{12} \Delta^{2} + \Delta c + d\right )$$



In [15]:

    
coeffs = solve([A-fm2, B-fm, C-f0, D-fp], [a,b,c,d], check=False)
coeffs









    Out[15]:





$$\left \{ a : \frac{1}{6 \Delta^{3}} \left(- 3 f_{i} + f_{{i+1}} + 3 f_{{i-1}} - f_{{i-2}}\right), \quad b : \frac{1}{2 \Delta^{2}} \left(- 2 f_{i} + f_{{i+1}} + f_{{i-1}}\right), \quad c : \frac{1}{24 \Delta} \left(15 f_{i} + 7 f_{{i+1}} - 27 f_{{i-1}} + 5 f_{{i-2}}\right), \quad d : \frac{13 f_{i}}{12} - \frac{f_{{i+1}}}{24} - \frac{f_{{i-1}}}{24}\right \}$$

and the pretty form of the polynomial



In [16]:

    
fc = f.subs(a,coeffs[a]).subs(b,coeffs[b]).subs(c,coeffs[c]).subs(d,coeffs[d])
fc









    Out[16]:





$$\frac{13 f_{i}}{12} - \frac{f_{{i+1}}}{24} - \frac{f_{{i-1}}}{24} + \frac{1}{24 \Delta} \left(x - x_{0}\right) \left(15 f_{i} + 7 f_{{i+1}} - 27 f_{{i-1}} + 5 f_{{i-2}}\right) + \frac{\left(x - x_{0}\right)^{2}}{2 \Delta^{2}} \left(- 2 f_{i} + f_{{i+1}} + f_{{i-1}}\right) + \frac{\left(x - x_{0}\right)^{3}}{6 \Delta^{3}} \left(- 3 f_{i} + f_{{i+1}} + 3 f_{{i-1}} - f_{{i-2}}\right)$$

this interpolant is symmetric about the $i-1/2$ interface---let's see the value there



In [18]:

    
fc.subs(x,x0-Rational(1,2)*dx)









    Out[18]:





$$\frac{7 f_{i}}{12} - \frac{f_{{i+1}}}{12} + \frac{7 f_{{i-1}}}{12} - \frac{f_{{i-2}}}{12}$$

Note that this is the interpolating polynomial used to find the interface states in PPM (Colella & Woodward 1984)

Quartic

Now we define a quartic polynomial that gives the correct averages over 5 zones, $x_{i-2}$, $x_{i-1}$, $x_i$, $x_{i+1}$, $x_{i+2}$, with zone averages $f_{i-2}$, $f_{i-1}$, $f_i$, $f_{i+1}$, $f_{i+2}$



In [18]:

    
a, b, c, d, e = symbols("a b c d e")
x0 = symbols("x0")
f = a*(x-x0)**4 + b*(x-x0)**3 + c*(x-x0)**2 + d*(x-x0) + e
f









    Out[18]:





$$a \left(x - x_{0}\right)^{4} + b \left(x - x_{0}\right)^{3} + c \left(x - x_{0}\right)^{2} + d \left(x - x_{0}\right) + e$$

Now we perform the integrals of $f(x)$ over each zone



In [24]:

    
fm2, fm, f0, fp, fp2 = symbols("f_{i-2} f_{i-1} f_i f_{i+1} f_{i+2}")
#xm32, xm12, xp12, xp32 = symbols("x_{i-3/2} x_{i-1/2} x_{i+1/2} x_{i+3/2}")
xm52 = x0 - Rational(5,2)*dx
xm32 = x0 - Rational(3,2)*dx
xm12 = x0 - Rational(1,2)*dx
xp12 = x0 + Rational(1,2)*dx
xp32 = x0 + Rational(3,2)*dx
xp52 = x0 + Rational(5,2)*dx

interfaces



In [25]:

    
xm52, xm32, xm12, xp12, xp32, xp52









    Out[25]:





$$\left ( - \frac{5 \Delta}{2} + x_{0}, \quad - \frac{3 \Delta}{2} + x_{0}, \quad - \frac{\Delta}{2} + x_{0}, \quad \frac{\Delta}{2} + x_{0}, \quad \frac{3 \Delta}{2} + x_{0}, \quad \frac{5 \Delta}{2} + x_{0}\right )$$



In [20]:

    
A = simplify(integrate(f/dx, (x, xm52, xm32)))
B = simplify(integrate(f/dx, (x, xm32, xm12)))
C = simplify(integrate(f/dx, (x, xm12, xp12)))
D = simplify(integrate(f/dx, (x, xp12, xp32)))
E = simplify(integrate(f/dx, (x, xp32, xp52)))

The analytic form of the constraints



In [21]:

    
A, B, C, D, E









    Out[21]:





$$\left ( \frac{1441 a}{80} \Delta^{4} - \frac{17 b}{2} \Delta^{3} + \frac{49 c}{12} \Delta^{2} - 2 \Delta d + e, \quad \frac{121 a}{80} \Delta^{4} - \frac{5 b}{4} \Delta^{3} + \frac{13 c}{12} \Delta^{2} - \Delta d + e, \quad \frac{\Delta^{4} a}{80} + \frac{\Delta^{2} c}{12} + e, \quad \frac{121 a}{80} \Delta^{4} + \frac{5 b}{4} \Delta^{3} + \frac{13 c}{12} \Delta^{2} + \Delta d + e, \quad \frac{1441 a}{80} \Delta^{4} + \frac{17 b}{2} \Delta^{3} + \frac{49 c}{12} \Delta^{2} + 2 \Delta d + e\right )$$

Now find the coefficients



In [22]:

    
coeffs = solve([A-fm2, B-fm, C-f0, D-fp, E-fp2], [a,b,c,d,e], check=False)
coeffs









    Out[22]:





$$\left \{ a : \frac{1}{24 \Delta^{4}} \left(6 f_{i} - 4 f_{{i+1}} + f_{{i+2}} - 4 f_{{i-1}} + f_{{i-2}}\right), \quad b : \frac{1}{12 \Delta^{3}} \left(- 2 f_{{i+1}} + f_{{i+2}} + 2 f_{{i-1}} - f_{{i-2}}\right), \quad c : \frac{1}{16 \Delta^{2}} \left(- 22 f_{i} + 12 f_{{i+1}} - f_{{i+2}} + 12 f_{{i-1}} - f_{{i-2}}\right), \quad d : \frac{1}{48 \Delta} \left(34 f_{{i+1}} - 5 f_{{i+2}} - 34 f_{{i-1}} + 5 f_{{i-2}}\right), \quad e : \frac{1067 f_{i}}{960} - \frac{29 f_{{i+1}}}{480} + \frac{3 f_{{i+2}}}{640} - \frac{29 f_{{i-1}}}{480} + \frac{3 f_{{i-2}}}{640}\right \}$$

and the pretty form of the polynomial



In [23]:

    
f.subs(a,coeffs[a]).subs(b,coeffs[b]).subs(c,coeffs[c]).subs(d,coeffs[d]).subs(e,coeffs[e])









    Out[23]:





$$\frac{1067 f_{i}}{960} - \frac{29 f_{{i+1}}}{480} + \frac{3 f_{{i+2}}}{640} - \frac{29 f_{{i-1}}}{480} + \frac{3 f_{{i-2}}}{640} + \frac{1}{48 \Delta} \left(x - x_{0}\right) \left(34 f_{{i+1}} - 5 f_{{i+2}} - 34 f_{{i-1}} + 5 f_{{i-2}}\right) + \frac{\left(x - x_{0}\right)^{2}}{16 \Delta^{2}} \left(- 22 f_{i} + 12 f_{{i+1}} - f_{{i+2}} + 12 f_{{i-1}} - f_{{i-2}}\right) + \frac{\left(x - x_{0}\right)^{3}}{12 \Delta^{3}} \left(- 2 f_{{i+1}} + f_{{i+2}} + 2 f_{{i-1}} - f_{{i-2}}\right) + \frac{\left(x - x_{0}\right)^{4}}{24 \Delta^{4}} \left(6 f_{i} - 4 f_{{i+1}} + f_{{i+2}} - 4 f_{{i-1}} + f_{{i-2}}\right)$$



In [ ]: