Finite Difference



In [2]:

    
import sympy
from sympy import init_printing, pprint, latex 
from sympy import binomial, Function, factorial, Symbol, symbols, Sum, S  
# S is for SingletonRegistry
from sympy import Rational as Rat
from sympy import Matrix, solve_linear_system

EY : 20160627 I have no idea what sympy.SingletonRegistry means despite the documentation for SingletonRegistry, which is sparse at best.



In [3]:

    
init_printing() # LaTeX printing



In [4]:

    
n = Symbol('n',integer=True,positive=True)
i = Symbol('i',integer=True,positive=True)
j = Symbol('j',integer=True,positive=True)
x,h,x_0 = symbols('x h x_0',real=True)
f = Function('f')



In [8]:

    
factorial(5)









    Out[8]:





120



In [18]:

    
C_n = Sum((-1)**i*binomial(n,i),(i,1,n))



In [24]:

    
(C_n.subs(n,3)).evalf()









    Out[24]:





-1.00000000000000

cf. From stackexchange How to do a symbolic taylor expansion of an unknown function $f(x)$ using sympy



In [78]:

    
f(x).series(x, x0=x_0, n=3)









    Out[78]:





$$f{\left (x_{0} \right )} + \left(x - x_{0}\right) \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \frac{1}{2} \left(x - x_{0}\right)^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}}  f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(\left(x - x_{0}\right)^{3}; x\rightarrowx_{0}\right)$$



In [32]:

    
pprint( f(x).series(x, x0=h, n=3) )









    



                                              ⎛  2        ⎞│                  
                                            2 ⎜ d         ⎟│                  
                                    (-h + x) ⋅⎜────(f(ξ₁))⎟│                  
                                              ⎜   2       ⎟│                  
                ⎛ d        ⎞│                 ⎝dξ₁        ⎠│ξ₁=h    ⎛        3
f(h) + (-h + x)⋅⎜───(f(ξ₁))⎟│     + ──────────────────────────── + O⎝(-h + x) 
                ⎝dξ₁       ⎠│ξ₁=h                2                            

        
        
        
        
       ⎞
; x → h⎠



In [33]:

    
latex(f(x).series(x, x0=h, n=3) )









    Out[33]:





'f{\\left (h \\right )} + \\left(- h + x\\right) \\left. \\frac{d}{d \\xi_{1}} f{\\left (\\xi_{1} \\right )} \\right|_{\\substack{ \\xi_{1}=h }} + \\frac{1}{2} \\left(- h + x\\right)^{2} \\left. \\frac{d^{2}}{d \\xi_{1}^{2}}  f{\\left (\\xi_{1} \\right )} \\right|_{\\substack{ \\xi_{1}=h }} + \\mathcal{O}\\left(\\left(- h + x\\right)^{3}; x\\rightarrowh\\right)'



In [37]:

    
FW = f(x+h).series(x+h, x0=x_0, n=3)
FW = FW.subs(x-x_0,0)



In [38]:

    
FW









    Out[38]:





$$f{\left (x_{0} \right )} + h \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \frac{h^{2}}{2} \left. \frac{d^{2}}{d \xi_{1}^{2}}  f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{3} + h^{2} x + h x^{2} + x^{3}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$



In [42]:

    
FW.args[0]









    Out[42]:





$$h \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }}$$



In [55]:

    
(f(x+h).series(x+h,x0=x_0,n=3)+f(x-h).series(x-h,x0=x_0,n=3)).subs(x-x_0,0)









    Out[55]:





$$2 f{\left (x_{0} \right )} + h^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}}  f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{3} + h^{2} x + h x^{2} + x^{3}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$



In [62]:

    
( f(x+h).series(x+h,x0=x_0,n=3) + f(x-h).series(x-h,x0=x_0,n=3) + \
f(x+2*h).series(x+2*h,x0=x_0,n=3) + f(x-2*h).series(x-2*h,x0=x_0,n=3) ).subs(x-x_0,0)









    Out[62]:





$$4 f{\left (x_{0} \right )} + 5 h^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}}  f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{3} + h^{2} x + h x^{2} + x^{3}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$



In [63]:

    
( f(x+h).series(x+h,x0=x_0,n=3) - f(x-h).series(x-h,x0=x_0,n=3) + \
f(x+2*h).series(x+2*h,x0=x_0,n=3) - f(x-2*h).series(x-2*h,x0=x_0,n=3) ).subs(x-x_0,0)









    Out[63]:





$$6 h \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{3} + h^{2} x + h x^{2} + x^{3}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$

finite-difference for the 1st order derivative ($n=1$), to the $p$th-order, $p=8$



In [77]:

    
p = 8
( Rat(4,5)*(f(x+h).series(x+h,x0=x_0,n=p) - f(x-h).series(x-h,x0=x_0,n=p)) + \
Rat(-1,5)*(f(x+2*h).series(x+2*h,x0=x_0,n=p)-f(x-2*h).series(x-2*h,x0=x_0,n=p)) + \
Rat(4,105)*(f(x+3*h).series(x+3*h,x0=x_0,n=p)-f(x-3*h).series(x-3*h,x0=x_0,n=p)) + \
Rat(-1,280)*(f(x+4*h).series(x+4*h,x0=x_0,n=p)-f(x-4*h).series(x-4*h,x0=x_0,n=p)) ).subs(x-x_0,0)









    Out[77]:





$$h \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{8} + h^{7} x + h^{6} x^{2} + h^{5} x^{3} + h^{4} x^{4} + h^{3} x^{5} + h^{2} x^{6} + h x^{7} + x^{8}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$



In [74]:

    
Rat(4,5)









    Out[74]:





$$\frac{4}{5}$$

centered difference approximation

$$-i_{\text{min}} = (d+p-1)/2 $$

$$i_{\text{max}} = \lceil (d+p-1)/2 \rceil $$

cf. Derivative Approximation by Finite Differences David Eberly



In [88]:

    
d = 2 
p = 8
sympy.floor(Rat(d+p-1,2)), Rat((d+p-1),2)









    Out[88]:





$$\left ( 4, \quad \frac{9}{2}\right )$$

First order derivatives $d=1$; let's aim for $-i_{\text{min}} = i_{\text{max}} = 2$ so we only need 2 neighbor values.



In [104]:

    
d = 1 
p = 5
sympy.floor(Rat(d+p-1,2)), Rat((d+p-1),2)









    Out[104]:





$$\left ( 2, \quad \frac{5}{2}\right )$$



In [106]:

    
d = 2 
p = 4
sympy.floor(Rat(d+p-1,2)), Rat((d+p-1),2)









    Out[106]:





$$\left ( 2, \quad \frac{5}{2}\right )$$



In [4]:

    
C_1,C_2,C_3,C_4 = symbols('C_1 C_2 C_3 C_4', integer=True)

From sympy documentation, Matrices (linear algebra), Basic Manipulation, I learned to construct matrices with the lambda function



In [117]:

    
Matrix(3,4, lambda i,j: (j+1)**(2*(i+1)+1))









    Out[117]:





$$\left[\begin{matrix}1 & 8 & 27 & 64\\1 & 32 & 243 & 1024\\1 & 128 & 2187 & 16384\end{matrix}\right]$$



In [116]:

    
Matrix(1,4,lambda i,j: 2*(j+1))









    Out[116]:





$$\left[\begin{matrix}2 & 4 & 6 & 8\end{matrix}\right]$$



In [138]:

    
C_nu_d1p6 = sympy.zeros(4,5) # d = 1, p = 6
C_nu_d1p6[:3,:4] = Matrix(3,4,lambda i,j: (j+1)**(2*(i+1)+1))
C_nu_d1p6[3,:4] = Matrix(1,4,lambda i,j: 2*(j+1))
C_nu_d1p6[3,4] = Rat(1)



In [139]:

    
C_nu_d1p6









    Out[139]:





$$\left[\begin{matrix}1 & 8 & 27 & 64 & 0\\1 & 32 & 243 & 1024 & 0\\1 & 128 & 2187 & 16384 & 0\\2 & 4 & 6 & 8 & 1\end{matrix}\right]$$



In [148]:

    
solve_linear_system( C_nu_d1p6, C_1,C_2,C_3,C_4 )









    Out[148]:





$$\left \{ C_{1} : \frac{4}{5}, \quad C_{2} : - \frac{1}{5}, \quad C_{3} : \frac{4}{105}, \quad C_{4} : - \frac{1}{280}\right \}$$

Also, another way is this:



In [142]:

    
C_nuM_d1p6 = sympy.zeros(4,4) # d = 1, p = 6
C_nuM_d1p6[:3,:4] = Matrix(3,4,lambda i,j: (j+1)**(2*(i+1)+1))
C_nuM_d1p6[3,:4] = Matrix(1,4,lambda i,j: 2*(j+1))
C_nuM_d1p6









    Out[142]:





$$\left[\begin{matrix}1 & 8 & 27 & 64\\1 & 32 & 243 & 1024\\1 & 128 & 2187 & 16384\\2 & 4 & 6 & 8\end{matrix}\right]$$



In [144]:

    
C_nuM_d1p6.inv() # and then read off the last column









    Out[144]:





$$\left[\begin{matrix}- \frac{61}{90} & \frac{29}{360} & - \frac{1}{360} & \frac{4}{5}\\\frac{169}{360} & - \frac{13}{180} & \frac{1}{360} & - \frac{1}{5}\\- \frac{1}{10} & \frac{1}{40} & - \frac{1}{840} & \frac{4}{105}\\\frac{7}{720} & - \frac{1}{360} & \frac{1}{5040} & - \frac{1}{280}\end{matrix}\right]$$



In [146]:

    
Matrix([[1,8],[2,4]]).inv() # d = 1, p=2









    Out[146]:





$$\left[\begin{matrix}- \frac{1}{3} & \frac{2}{3}\\\frac{1}{6} & - \frac{1}{12}\end{matrix}\right]$$



In [4]:

    
C_nu_d1p4 = sympy.zeros(3,4) # d = 1, p = 4
C_nu_d1p4[:2,:3] = Matrix(2,3,lambda i,j: (j+1)**(2*(i+1)+1))
C_nu_d1p4[2,:3] = Matrix(1,3,lambda i,j: 2*(j+1))
C_nu_d1p4[2,3] = Rat(1)



In [5]:

    
C_nu_d1p4









    Out[5]:





$$\left[\begin{matrix}1 & 8 & 27 & 0\\1 & 32 & 243 & 0\\2 & 4 & 6 & 1\end{matrix}\right]$$



In [8]:

    
solve_linear_system( C_nu_d1p4, C_1,C_2,C_3 )









    Out[8]:





$$\left \{ C_{1} : \frac{3}{4}, \quad C_{2} : - \frac{3}{20}, \quad C_{3} : \frac{1}{60}\right \}$$

Indeed,



In [9]:

    
p = 6
( Rat(3,4)*(f(x+h).series(x+h,x0=x_0,n=p) - f(x-h).series(x-h,x0=x_0,n=p)) + \
Rat(-3,20)*(f(x+2*h).series(x+2*h,x0=x_0,n=p)-f(x-2*h).series(x-2*h,x0=x_0,n=p)) + \
Rat(1,60)*(f(x+3*h).series(x+3*h,x0=x_0,n=p)-f(x-3*h).series(x-3*h,x0=x_0,n=p))).subs(x-x_0,0)









    Out[9]:





$$h \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{6} + h^{5} x + h^{4} x^{2} + h^{3} x^{3} + h^{2} x^{4} + h x^{5} + x^{6}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$

$$ \boxed{ h \frac{d}{dx}f(x) = \frac{3}{4} ( f(x+h) - f(x-h) ) + \frac{-3}{20}( f(x+2h) - f(x-2h)) + \frac{1}{60} (f(x+3h)-f(x-3h)) +O(h^6) = \sum_{\nu = 1}^3 C_{\nu} (f(x+ \nu h) - f(x - \nu h) ) + O(h^6) } $$

with

$C_1 = \frac{3}{4}$
$C_2 = \frac{-3}{20}$
$C_3 = \frac{1}{60}$

Let's generalize this procedure:



In [60]:

    
def make_C_nu_dp(NU):
    C_nu_dp = sympy.zeros(NU,NU+1)
    C_nu_dp[:NU-1,:NU] = Matrix(NU-1,NU,lambda i,j: (j+1)**(2*(i+1)+1))
    #C_nu_dp[NU-1,:NU]  = Matrix(1,NU,lambda i,j:2*(j+1))
    C_nu_dp[NU-1,:NU]  = Matrix(1,NU,lambda i,j:(j+1))
    C_nu_dp[NU-1,NU]   = Rat(1,2)

    # sanity check
    pprint( C_nu_dp)
    Cnus = symbols(' '.join( ['C_'+str(nu)for nu in range(1,NU+1)]), real=True)
    return solve_linear_system( C_nu_dp, *Cnus )



In [61]:

    
make_C_nu_dp(2)









    



⎡1  8   0 ⎤
⎢         ⎥
⎣1  2  1/2⎦






    Out[61]:





{C_1: 2/3, C_2: -1/12}



In [70]:

    
# sanity check
p=2
(Rat(2,3)*(f(x+h).series(x+h,x0=x_0,n=p) - f(x-h).series(x-h,x0=x_0,n=p)) + \
Rat(-1,12)*(f(x+2*h).series(x+2*h,x0=x_0,n=p)-f(x-2*h).series(x-2*h,x0=x_0,n=p))).simplify()









    Out[70]:





$$h \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{2} + h x + x^{2}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$



In [71]:

    
make_C_nu_dp(3)









    



⎡1  8   27    0 ⎤
⎢               ⎥
⎢1  32  243   0 ⎥
⎢               ⎥
⎣1  2    3   1/2⎦






    Out[71]:





$$\left \{ C_{1} : \frac{3}{4}, \quad C_{2} : - \frac{3}{20}, \quad C_{3} : \frac{1}{60}\right \}$$



In [72]:

    
make_C_nu_dp(4)









    



⎡1   8    27    64     0 ⎤
⎢                        ⎥
⎢1  32   243   1024    0 ⎥
⎢                        ⎥
⎢1  128  2187  16384   0 ⎥
⎢                        ⎥
⎣1   2    3      4    1/2⎦






    Out[72]:





$$\left \{ C_{1} : \frac{4}{5}, \quad C_{2} : - \frac{1}{5}, \quad C_{3} : \frac{4}{105}, \quad C_{4} : - \frac{1}{280}\right \}$$

Let's attempt to further generalize this procedure for $d=1$.

Consider the following easy examples for $p=1$.

Now $$ f(x \pm h) = f(x) + (\frac{ \partial f}{ \partial x}(x) )( \pm h) + \mathcal{O}(h^2) $$
$$ f(x+h) - f(x-h) = 2h \frac{ \partial f}{ \partial x} + \mathcal{O}(h^2) \text{ or } \frac{ \partial f}{ \partial x} = \frac{ f(x+h) - f(x-h) }{2h} + \mathcal{O}(h) $$

Likewise, for so-called forward and backward difference operators: $$ \frac{f(x+h) - f(x)}{h} = \frac{ \partial f}{ \partial x} $$ $$ \frac{f(x) - f(x-h)}{h} = \frac{ \partial f}{ \partial x} $$ Conjecture that $$ \boxed{ \frac{ \partial f}{ \partial x} = \sum^N_{\nu = 1} C_{\nu} (f(x+ \nu h) - f(x-\nu h) ) \left( \frac{1}{h} \right) } $$

and so doing the Taylor expansion, $$ f(x \pm \nu h) =\sum_{j=0}^{p} \frac{ f^{(j)}(x) }{j!} ( \pm \nu h)^j + \mathcal{O}(h^{p+1}) = \sum_{j'=0}^{p'} \frac{ f^{(2j')}(x) }{ (2j')! } (\nu h)^{2j'} + \sum_{j'=0}^{p'} \frac{ f^{(2j'+1) }(x) }{ (2j'+1)!} (\nu h)^{2j'+1} (\pm 1) + \mathcal{O}(h^{p+1} ) $$ where $p := 2p'+1$.

Clearly, for $f(x+\nu h) - f(x-\nu h)$, all the "even" terms in the Taylor series expansion cancel each other out, leaving $$ f(x+\nu h) - f(x-\nu h) = \sum_{j'=0}^{p'} \frac{ 2f^{(2j'+1)}(x) }{(2j'+1)!} (\nu h)^{2j'+1} + O(h^{p+1}) = 2\nu h \frac{ \partial f}{\partial x} + \sum_{j'=1}^{p'} \frac{2(\nu h)^{2j' +1} }{ (2j'+1)! } f^{(2j'+1) }(x) + \mathcal{O}(h^{p+1} ) $$

and so with our previous assumption (caveat),

$$ \sum_{\nu =1}^{N} C_{\nu} ( f(x+ \nu h) - f(x-\nu h)) \frac{1}{h} = \sum_{\nu =1}^N C_{\nu} \left( 2\nu \frac{ \partial f}{ \partial x} + \sum_{j'=1}^{p'} \frac{ 2\nu ^{2j'+1} }{ (2j'+1)! } f^{(2j+1) }(x) h^{2j'} \right) + \mathcal{O}(h^{p} ) $$

Then we obtain a system of linear equations that we'd want to solve:

$$ \sum_{\nu =1}^N 2C_{\nu} \nu = 1 $$

$$ \sum_{j'=1}^{p'} \left( \sum_{\nu =1}^N C_{\nu} \frac{2\nu^{2j'+1} }{ (2j'+1)! } \right) f^{(2j+1) }(x) h^{2j'} = 0 \text{ or } \left( \sum_{\nu =1}^N C_{\nu} \frac{ 2 \nu^{2j'+1} }{ (2j'+1)! } \right) = 0 \qquad \, \forall \, j' = 1 \dots p' $$

Thus, we confirm the make_C_nu_dp(NU) Python function for $d=1$.

Backward difference operator, in the case of $d=1$, first-order derivative

For the backward difference operator, then consider the following:

$$ \frac{ \partial f}{ \partial x} = \sum_{ \nu =0}^N \frac{1}{h} C_{\nu} f(x-\nu h) = \sum_{\nu = 0}^N \frac{C_{\nu} }{ h} \left[ f(x) + \left( \frac{ \partial f}{ \partial x} \right) (-\nu h) + \sum_{j=2}^p \frac{ f^{(j)}(x) }{j!} (-\nu h)^j \right] + \mathcal{O}(h^p) $$

which leads to a system of linear equations:

$$ \frac{1}{h} \sum_{\nu =0}^N C_{\nu} = 0 $$$$ - \sum_{\nu =0}^N C_{\nu } \nu = 1 $$$$ \sum_{j=2}^p \frac{(-h)^jf^{(j)}(x)}{ hj!} \sum_{ \nu =0}^N C_{\nu} \nu^j = 0 $$

and I didn't make it explicit before, and neither here, but you can divide both sides of the equation by constant values that don't depend on index $j$ (or $\nu$) and obtain the system of linear equations to solve.



In [126]:

    
NU = 3
C_nueqns = sympy.zeros(NU+1,NU+2)
C_nueqns[0,:NU+1] = Matrix(1,NU+1,lambda i,j: 1)
C_nueqns[1,:NU+1] = Matrix(1,NU+1,lambda i,j:j)
C_nueqns[1,NU+1] = -1
if NU > 1:
    C_nueqns[2:NU+1,:NU+1] = Matrix(NU-1,NU+1,lambda i,j: j**(i+2) )
pprint(C_nueqns)
Cnus = symbols(' '.join( ['C_'+str(nu)for nu in range(NU+1)]), real=True)
pprint(Cnus)
solve_linear_system( C_nueqns,*Cnus)









    



⎡1  1  1  1   0 ⎤
⎢               ⎥
⎢0  1  2  3   -1⎥
⎢               ⎥
⎢0  1  4  9   0 ⎥
⎢               ⎥
⎣0  1  8  27  0 ⎦
(C₀, C₁, C₂, C₃)






    Out[126]:





$$\left \{ C_{0} : \frac{11}{6}, \quad C_{1} : -3, \quad C_{2} : \frac{3}{2}, \quad C_{3} : - \frac{1}{3}\right \}$$



In [127]:

    
def make_backwardC_nu(NU):
    C_nueqns = sympy.zeros(NU+1,NU+2)
    C_nueqns[0,:NU+1]  = Matrix(1,NU+1,lambda i,j: 1)

    C_nueqns[1,:NU+1]  = Matrix(1,NU+1,lambda i,j: j)
    C_nueqns[1,NU+1]   = -1

    if NU > 1:
        C_nueqns[2:NU+1,:NU+1] = Matrix(NU-1,NU+1,lambda i,j: (j)**(i+2) )

    # sanity check
    pprint( C_nueqns)    
        
    Cnus = symbols(' '.join( ['C_'+str(nu)for nu in range(NU+1)]), real=True)
    return solve_linear_system( C_nueqns, *Cnus )



In [122]:

    
make_backwardC_nu(1)









    



⎡1  1  0 ⎤
⎢        ⎥
⎣0  1  -1⎦






    Out[122]:





$$\left \{ C_{0} : 1, \quad C_{1} : -1\right \}$$



In [128]:

    
make_backwardC_nu(2)









    



⎡1  1  1  0 ⎤
⎢           ⎥
⎢0  1  2  -1⎥
⎢           ⎥
⎣0  1  4  0 ⎦






    Out[128]:





$$\left \{ C_{0} : \frac{3}{2}, \quad C_{1} : -2, \quad C_{2} : \frac{1}{2}\right \}$$



In [168]:

    
# sanity check
f(x).series(x,x0=x_0,n=4).subs(x-x_0,0)
f(x-h).series(x-h,x0=x_0,n=4).subs(x-x_0,0)
f(x-2*h).series(x-2*h,x0=x_0,n=4).subs(x-x_0,0)

p=4
( Rat(3,2)*f(x).series(x,x0=x_0,n=p).subs(x-x_0,0) + \
  Rat(-2)*( f(x-h).series(x-h,x0=x_0,n=p) ).subs(x-x_0,0) + \
Rat(1,2)*( f(x-2*h).series(x-2*h,x0=x_0,n=p)).subs(x-x_0,0) )









    Out[168]:





$$h \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} - \frac{h^{3}}{3} \left. \frac{d^{3}}{d \xi_{1}^{3}}  f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{4} + h^{3} x + h^{2} x^{2} + h x^{3} + x^{4}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$



In [137]:

    
make_backwardC_nu(3)









    



⎡1  1  1  1   0 ⎤
⎢               ⎥
⎢0  1  2  3   -1⎥
⎢               ⎥
⎢0  1  4  9   0 ⎥
⎢               ⎥
⎣0  1  8  27  0 ⎦






    Out[137]:





$$\left \{ C_{0} : \frac{11}{6}, \quad C_{1} : -3, \quad C_{2} : \frac{3}{2}, \quad C_{3} : - \frac{1}{3}\right \}$$



In [169]:

    
# sanity check
p=5
( Rat(11,6)*f(x).series(x,x0=x_0,n=p).subs(x-x_0,0) + \
  Rat(-3)*( f(x-h).series(x-h,x0=x_0,n=p) ).subs(x-x_0,0) + \
Rat(3,2)*( f(x-2*h).series(x-2*h,x0=x_0,n=p)).subs(x-x_0,0) + \
Rat(-1,3)*( f(x-3*h).series(x-3*h,x0=x_0,n=p)).subs(x-x_0,0) 
)









    Out[169]:





$$h \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} - \frac{h^{4}}{4} \left. \frac{d^{4}}{d \xi_{1}^{4}}  f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{5} + h^{4} x + h^{3} x^{2} + h^{2} x^{3} + h x^{4} + x^{5}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$



In [138]:

    
make_backwardC_nu(4)









    



⎡1  1  1   1    1   0 ⎤
⎢                     ⎥
⎢0  1  2   3    4   -1⎥
⎢                     ⎥
⎢0  1  4   9   16   0 ⎥
⎢                     ⎥
⎢0  1  8   27  64   0 ⎥
⎢                     ⎥
⎣0  1  16  81  256  0 ⎦






    Out[138]:





$$\left \{ C_{0} : \frac{25}{12}, \quad C_{1} : -4, \quad C_{2} : 3, \quad C_{3} : - \frac{4}{3}, \quad C_{4} : \frac{1}{4}\right \}$$



In [171]:

    
# sanity check
p=6
( Rat(25,12)*f(x).series(x,x0=x_0,n=p).subs(x-x_0,0) + \
  Rat(-4)*( f(x-h).series(x-h,x0=x_0,n=p) ).subs(x-x_0,0) + \
Rat(3)*( f(x-2*h).series(x-2*h,x0=x_0,n=p)).subs(x-x_0,0) + \
Rat(-4,3)*( f(x-3*h).series(x-3*h,x0=x_0,n=p)).subs(x-x_0,0) + \
Rat(1,4)*( f(x-4*h).series(x-4*h,x0=x_0,n=p)).subs(x-x_0,0) 
)









    Out[171]:





$$h \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} - \frac{h^{5}}{5} \left. \frac{d^{5}}{d \xi_{1}^{5}}  f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{6} + h^{5} x + h^{4} x^{2} + h^{3} x^{3} + h^{2} x^{4} + h x^{5} + x^{6}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$

For $d=2$



In [164]:

    
C_nu_d2p4 = sympy.zeros(2,3) # d = 2, p = 4
C_nu_d2p4[:1,:2] = Matrix(1,2,lambda i,j: (j+1)**(2*(i+2)))
C_nu_d2p4[1,:2] = Matrix(1,2,lambda i,j: 2*(j+1)**2)
C_nu_d2p4[1,2] = Rat(1)



In [165]:

    
C_nu_d2p4









    Out[165]:





$$\left[\begin{matrix}1 & 16 & 0\\2 & 8 & 1\end{matrix}\right]$$



In [166]:

    
solve_linear_system( C_nu_d2p4, C_1,C_2)









    Out[166]:





$$\left \{ C_{1} : \frac{2}{3}, \quad C_{2} : - \frac{1}{24}\right \}$$



In [168]:

    
C_nu_d2p6 = sympy.zeros(3,4) # d = 2, p = 6
C_nu_d2p6[:2,:3] = Matrix(2,3,lambda i,j: (j+1)**(2*(i+2)))
C_nu_d2p6[2,:3] = Matrix(1,3,lambda i,j: 2*(j+1)**2)
C_nu_d2p6[2,3] = Rat(1)



In [169]:

    
C_nu_d2p6









    Out[169]:





$$\left[\begin{matrix}1 & 16 & 81 & 0\\1 & 64 & 729 & 0\\2 & 8 & 18 & 1\end{matrix}\right]$$



In [170]:

    
solve_linear_system( C_nu_d2p6, C_1,C_2,C_3)









    Out[170]:





$$\left \{ C_{1} : \frac{3}{4}, \quad C_{2} : - \frac{3}{40}, \quad C_{3} : \frac{1}{180}\right \}$$

Indeed,



In [172]:

    
p = 4
( Rat(2,3)*(f(x+h).series(x+h,x0=x_0,n=p) + f(x-h).series(x-h,x0=x_0,n=p) -2*f(x_0)) + \
Rat(-1,24)*(f(x+2*h).series(x+2*h,x0=x_0,n=p)+f(x-2*h).series(x-2*h,x0=x_0,n=p) -2*f(x_0))).subs(x-x_0,0)









    Out[172]:





$$\frac{h^{2}}{2} \left. \frac{d^{2}}{d \xi_{1}^{2}}  f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{4} + h^{3} x + h^{2} x^{2} + h x^{3} + x^{4}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$



In [175]:

    
p = 6
( Rat(3,4)*(f(x+h).series(x+h,x0=x_0,n=p) + f(x-h).series(x-h,x0=x_0,n=p)-2*f(x_0)) + \
Rat(-3,40)*(f(x+2*h).series(x+2*h,x0=x_0,n=p)+f(x-2*h).series(x-2*h,x0=x_0,n=p)-2*f(x_0)) + \
Rat(1,180)*(f(x+3*h).series(x+3*h,x0=x_0,n=p)+f(x-3*h).series(x-3*h,x0=x_0,n=p)-2*f(x_0))).subs(x-x_0,0)









    Out[175]:





$$\frac{h^{2}}{2} \left. \frac{d^{2}}{d \xi_{1}^{2}}  f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{6} + h^{5} x + h^{4} x^{2} + h^{3} x^{3} + h^{2} x^{4} + h x^{5} + x^{6}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$

Central difference operator, in the case of $d=2$, 2nd-order derivative

Let's attempt to further generalize this procedure for $d=2$.

Consider the following easy examples for $p=1$.

Now $$ f(x \pm h) = f(x) + (\frac{ \partial f}{ \partial x}(x) )( \pm h) + ( \frac{1}{2} \frac{ \partial^2 f}{ \partial x^2}(x) )h^2 + \mathcal{O}(h^2) $$
$$ f(x+h) + f(x-h) = 2f(x) + h^2 \frac{ \partial^2 f}{ \partial x^2} + \mathcal{O}(h^3) \text{ or } \frac{ \partial^2 f}{ \partial x^2} = \frac{ f(x+h) + f(x-h) - 2f(x) }{h^2} + \mathcal{O}(h) $$

And so

$$ f(x \pm \nu h) = f(x) + (\frac{ \partial f}{ \partial x}(x) )( \pm \nu h) + (\frac{1}{2} \frac{ \partial^2 f}{ \partial x^2}(x) ) \nu^2 h^2 + \mathcal{O}(h^2) $$

$$ f(x+\nu h) + f(x-\nu h) = 2f(x) + \nu^2 h^2 \frac{ \partial^2 f}{ \partial x^2} + \mathcal{O}(h^3) \text{ or } \frac{ \partial^2 f}{ \partial x^2} = \frac{ f(x+\nu h) + f(x-\nu h) - 2f(x) }{\nu^2 h^2} + \mathcal{O}(h) $$

Conjecture that $$ \boxed{ \frac{ \partial^2 f}{ \partial x^2} = \sum^N_{\nu = 1} C_{\nu} (f(x+ \nu h) + f(x-\nu h) -2f(x) ) \left( \frac{1}{h^2} \right) } $$ (yes, the $\frac{1}{\nu^2}$ factor was swallowed into $C_{\nu}$.)

Clearly, for $f(x+\nu h) + f(x-\nu h)$, all the "odd" terms in the Taylor series expansion cancel each other out, leaving $$ f(x+\nu h) + f(x-\nu h) -2f(x)= \sum_{j'=1}^{p'} \frac{ 2f^{(2j')}(x) }{(2j')!} (\nu h)^{2j'} + O(h^{p+1}) = \frac{ 2f''(x)}{2} (\nu h)^2 + \sum_{j'=2}^{p'} \frac{2(\nu h)^{2j' } }{ (2j')! } f^{(2j') }(x) + \mathcal{O}(h^{p+1} ) $$

and so with our previous assumption (caveat),

$$ \sum_{\nu =1}^{N} C_{\nu} ( f(x+ \nu h) + f(x-\nu h) -2f(x)) \frac{1}{h^2} = \sum_{\nu =1}^N C_{\nu} \left( \nu^2 f''(x)h^2 + \sum_{j'=2}^{p'} \frac{ 2\nu ^{2j'} }{ (2j')! } f^{(2j') }(x) h^{2j'} \right) \frac{1}{h^2} + \mathcal{O}(h^{p} ) = \sum_{\nu =1}^N \left[ C_{\nu} \nu^2 f''(x) + C_{\nu} \sum_{j'=2}^{p'} \frac{ 2\nu^{2j'} }{ (2j')!} f^{(2j') }(x)h^{2(j'-1) } \right] + \mathcal{O}(h^p) $$

Then we obtain a system of linear equations that we'd want to solve:

$$ \sum_{\nu =1}^N \nu^2 C_{\nu} = 1 $$

$$ \sum_{j'=2}^{p'} \left( \sum_{\nu =1}^N C_{\nu} \frac{ 2\nu^{2j'} }{ (2j')! } \right)= 0 $$



In [36]:

    
def make_C_nu_d2p(NU):
    C_nu_dp = sympy.zeros(NU,NU+1)
    C_nu_dp[:NU-1,:NU] = Matrix(NU-1,NU,lambda i,j: (j+1)**(2*(i+2)))

    C_nu_dp[NU-1,:NU]  = Matrix(1,NU,lambda i,j:(j+1)**2)
    C_nu_dp[NU-1,NU]   = Rat(1)

    # sanity check
    pprint( C_nu_dp)
    Cnus = symbols(' '.join( ['C_'+str(nu)for nu in range(1,NU+1)]), real=True)
    return solve_linear_system( C_nu_dp, *Cnus )



In [37]:

    
make_C_nu_d2p(2)









    



⎡1  16  0⎤
⎢        ⎥
⎣1  4   1⎦






    Out[37]:





$$\left \{ C_{1} : \frac{4}{3}, \quad C_{2} : - \frac{1}{12}\right \}$$



In [38]:

    
make_C_nu_d2p(3)









    



⎡1  16  81   0⎤
⎢             ⎥
⎢1  64  729  0⎥
⎢             ⎥
⎣1  4    9   1⎦






    Out[38]:





$$\left \{ C_{1} : \frac{3}{2}, \quad C_{2} : - \frac{3}{20}, \quad C_{3} : \frac{1}{90}\right \}$$



In [39]:

    
make_C_nu_d2p(4)









    



⎡1  16    81    256   0⎤
⎢                      ⎥
⎢1  64   729   4096   0⎥
⎢                      ⎥
⎢1  256  6561  65536  0⎥
⎢                      ⎥
⎣1   4    9     16    1⎦






    Out[39]:





$$\left \{ C_{1} : \frac{8}{5}, \quad C_{2} : - \frac{1}{5}, \quad C_{3} : \frac{8}{315}, \quad C_{4} : - \frac{1}{560}\right \}$$

Indeed



In [40]:

    
p = 6
( Rat(4,3)*(f(x+h).series(x+h,x0=x_0,n=p) + f(x-h).series(x-h,x0=x_0,n=p) -2*f(x_0)) + \
Rat(-1,12)*(f(x+2*h).series(x+2*h,x0=x_0,n=p)+f(x-2*h).series(x-2*h,x0=x_0,n=p) -2*f(x_0))).subs(x-x_0,0)









    Out[40]:





$$h^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}}  f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{6} + h^{5} x + h^{4} x^{2} + h^{3} x^{3} + h^{2} x^{4} + h x^{5} + x^{6}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$



In [41]:

    
p = 8
( Rat(3,2)*(f(x+h).series(x+h,x0=x_0,n=p) + f(x-h).series(x-h,x0=x_0,n=p) -2*f(x_0)) + \
Rat(-3,20)*(f(x+2*h).series(x+2*h,x0=x_0,n=p)+f(x-2*h).series(x-2*h,x0=x_0,n=p) -2*f(x_0)) + \
Rat(1,90)*(f(x+3*h).series(x+3*h,x0=x_0,n=p)+f(x-3*h).series(x-3*h,x0=x_0,n=p) -2*f(x_0)) ).subs(x-x_0,0)









    Out[41]:





$$h^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}}  f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{8} + h^{7} x + h^{6} x^{2} + h^{5} x^{3} + h^{4} x^{4} + h^{3} x^{5} + h^{2} x^{6} + h x^{7} + x^{8}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$



In [42]:

    
p = 10
( Rat(8,5)*(f(x+h).series(x+h,x0=x_0,n=p) + f(x-h).series(x-h,x0=x_0,n=p) -2*f(x_0)) + \
Rat(-1,5)*(f(x+2*h).series(x+2*h,x0=x_0,n=p)+f(x-2*h).series(x-2*h,x0=x_0,n=p) -2*f(x_0)) + \
Rat(8,315)*(f(x+3*h).series(x+3*h,x0=x_0,n=p)+f(x-3*h).series(x-3*h,x0=x_0,n=p) -2*f(x_0)) + \
 Rat(-1,560)*(f(x+4*h).series(x+4*h,x0=x_0,n=p)+f(x-4*h).series(x-4*h,x0=x_0,n=p) -2*f(x_0))
).subs(x-x_0,0)









    Out[42]:





$$h^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}}  f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{10} + h^{9} x + h^{8} x^{2} + h^{7} x^{3} + h^{6} x^{4} + h^{5} x^{5} + h^{4} x^{6} + h^{3} x^{7} + h^{2} x^{8} + h x^{9} + x^{10}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$

finite difference sanity check on Python numpy



In [5]:

    
import numpy as np



In [13]:

    
import matplotlib
import pylab
import numpy as np



In [14]:

    
%matplotlib inline
%config InlineBackend.figure_formats = {'svg','png'}



In [28]:

    
l_x = 2.* np.pi
l_y = 2.* np.pi
l_z = 2.* np.pi
L_x = 64 
L_y = 64 
L_z = 64
f_grid3d = np.zeros((L_x,L_y,L_z))
df_grid3d_result = np.zeros((L_x,L_y,L_z))
h_x = l_x/L_x; h_y = l_y/L_y; h_z = l_z/L_z



In [33]:

    
for i_z in range(L_z):
    for i_y in range(L_y):
        for i_x in range(L_x):
            f_grid3d[i_x][i_y][i_z] = np.sin( i_x*h_x)



In [26]:

    
f_grid3d[:,0,0].shape









    Out[26]:





(64,)



In [27]:

    
pylab.plot( f_grid3d[:,5,5] )









    Out[27]:





[<matplotlib.lines.Line2D at 0x7fe5ad8a5590>]



In [37]:

    
RAD = 3



In [38]:

    
for i_z in range(L_z-RAD):
    for i_y in range(L_y-RAD):
        for i_x in range(L_x-RAD):
#            if (i_z + RAD >= L_z): 
            df_grid3d_result[i_x][i_y][i_z] = \
                (3./(4.*h_x))*(f_grid3d[i_x+1][i_y][i_z]-f_grid3d[i_x-1][i_y][i_z]) + \
                (-3./(20.*h_x))*(f_grid3d[i_x+2][i_y][i_z]-f_grid3d[i_x-2][i_y][i_z]) + \
                (1./(60.*h_x))*(f_grid3d[i_x+3][i_y][i_z]-f_grid3d[i_x-3][i_y][i_z])



In [40]:

    
pylab.plot( df_grid3d_result[:,5,5] )









    Out[40]:





[<matplotlib.lines.Line2D at 0x7fe5ad7de150>]



In [42]:

    
pylab.plot( df_grid3d_result[:,0,0] )









    Out[42]:





[<matplotlib.lines.Line2D at 0x7fe5ad6c9610>]



In [35]:

    
f_grid3d[-2][0][0]









    Out[35]:





-0.19509032201612872



In [ ]:

    
p = 6
( Rat(3,4)*(f(x+h).series(x+h,x0=x_0,n=p) - f(x-h).series(x-h,x0=x_0,n=p)) + \
Rat(-3,20)*(f(x+2*h).series(x+2*h,x0=x_0,n=p)-f(x-2*h).series(x-2*h,x0=x_0,n=p)) + \
Rat(1,60)*(f(x+3*h).series(x+3*h,x0=x_0,n=p)-f(x-3*h).series(x-3*h,x0=x_0,n=p))).subs(x-x_0,0)

Instantons

cf. Valery Rubakov. Classical Theory of Gauge Fields. Princeton, 2002. Ch. 7. 7.1 Kink



In [182]:

    
t = Symbol('t',real=True)
phi1pl1 = Symbol(r'\varphi')# phi1pl1 denote phi 1+1 spacetime dimension; pl for "plus"
phi1pl1 = Function(r'\varphi')(t,x)



In [183]:

    
phi1pl1









    Out[183]:





$$\varphi{\left (t,x \right )}$$



In [191]:

    
phi1pl1.series(x=x,x0=x_0,n=3)









    Out[191]:





$$\varphi{\left (t,x_{0} \right )} + \left(x - x_{0}\right) \left. \frac{\partial}{\partial \xi_{2}} \varphi{\left (t,\xi_{2} \right )} \right|_{\substack{ \xi_{2}=x_{0} }} + \frac{1}{2} \left(x - x_{0}\right)^{2} \left. \frac{\partial^{2}}{\partial \xi_{2}^{2}}  \varphi{\left (t,\xi_{2} \right )} \right|_{\substack{ \xi_{2}=x_{0} }} + \mathcal{O}\left(\left(x - x_{0}\right)^{3}; x\rightarrowx_{0}\right)$$



In [193]:

    
((phi1pl1.diff(t)*phi1pl1.diff(t)) + (phi1pl1.diff(x)*phi1pl1.diff(x))).series(x=x,x0=x_0,n=3)









    Out[193]:





$$\left. \frac{\partial}{\partial x} \varphi{\left (t,x \right )} \right|_{\substack{ x=x + x_{0} }}^{2} + \frac{\partial}{\partial t} \varphi{\left (t,x_{0} \right )}^{2} + 2 \frac{\partial}{\partial t}\left(\frac{1}{2} \left(x - x_{0}\right)^{2} \left. \frac{\partial^{2}}{\partial \xi_{2}^{2}}  \varphi{\left (t,\xi_{2} \right )} \right|_{\substack{ \xi_{2}=x_{0} }}\right) \frac{\partial}{\partial t} \varphi{\left (t,x_{0} \right )} + 2 \frac{\partial}{\partial t}\left(\left(x - x_{0}\right) \left. \frac{\partial}{\partial \xi_{2}} \varphi{\left (t,\xi_{2} \right )} \right|_{\substack{ \xi_{2}=x_{0} }}\right) \frac{\partial}{\partial t} \varphi{\left (t,x_{0} \right )} + \frac{\partial}{\partial t}\left(\left(x - x_{0}\right) \left. \frac{\partial}{\partial \xi_{2}} \varphi{\left (t,\xi_{2} \right )} \right|_{\substack{ \xi_{2}=x_{0} }}\right)^{2} + \mathcal{O}\left(\left(x - x_{0}\right)^{3}; x\rightarrowx_{0}\right)$$



In [199]:

    
lambdavariablesymbol = Symbol(r'\lambda')
v = Symbol('v')
mu = Symbol(r'\mu')
m = Symbol('m')
lambdavariablesymbol/Rat(4)*(phi1pl1*phi1pl1-v**2)**2









    Out[199]:





$$\frac{\lambda}{4} \left(- v^{2} + \varphi^{2}{\left (t,x \right )}\right)^{2}$$



In [201]:

    
(lambdavariablesymbol/Rat(4)*(phi1pl1*phi1pl1-v**2)**2).series(x=x,x0=x_0,n=2).subs(v, \
                                                        mu/sympy.sqrt( lambdavariablesymbol))









    Out[201]:





$$\frac{\mu^{4}}{4 \lambda} + \left(x - x_{0}\right) \left(\lambda \varphi^{3}{\left (t,x_{0} \right )} \left. \frac{\partial}{\partial \xi_{2}} \varphi{\left (t,\xi_{2} \right )} \right|_{\substack{ \xi_{2}=x_{0} }} - \mu^{2} \varphi{\left (t,x_{0} \right )} \left. \frac{\partial}{\partial \xi_{2}} \varphi{\left (t,\xi_{2} \right )} \right|_{\substack{ \xi_{2}=x_{0} }}\right) - \frac{\mu^{2}}{2} \varphi^{2}{\left (t,x_{0} \right )} + \frac{\lambda}{4} \varphi^{4}{\left (t,x_{0} \right )} + \mathcal{O}\left(\left(x - x_{0}\right)^{2}; x\rightarrowx_{0}\right)$$

Dealing with Boundary conditions (with Finite Difference)

I found that I needed to represent boundary conditions that are "discretized" i.e. "numerical boundary condition" - for example on pp. 148 of Cebeci, Shao, Kafyeke, Laurendeau, (E5.2.3), I need to figure out if this is a derivative or second-order derivative:

$u_I^n = 2u^n_{I-1} - u^n_{I-2}$



In [10]:

    
(f(x).series(x,x0=x_0,n=3)- 2 * f(x-h).series(x-h,x0=x_0,n=3) + \
 f(x-2*h).series(x-2*h,x0=x_0,n=3)).subs(x-x_0,0)









    Out[10]:





$$h^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}}  f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{3} + h^{2} x + h x^{2} + x^{3}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$

So indeed this is essentially

$ \frac{ \partial^2 u}{ \partial x^2} = 0$

Also consider pp. 338 of Cebeci, Shao, Kafyeke, Laurendeau (2005), namely (11.5.7a), (11.5.7b), which purports the second-order accuracy of a boundary condition, involving a first-order partial derivative:

$ \frac{ \partial u^x }{ \partial y } = 0 $ and

$\frac{\partial p }{ \partial y} = 0 $

cf. on the "symmetry plane" of $y=H$ (these boundary conditions, $\lbrace (x,y) \in \mathbb{R}^2 | y=H \rbrace$ is a submanifold of $\mathbb{R}^2$ as it's a level set or from the locally imbedded theorem (from differential topology))



In [6]:

    
p=3
( Rat(3,2)*f(x).series(x,x0=x_0,n=p).subs(x-x_0,0) + \
  Rat(-2)*( f(x-h).series(x-h,x0=x_0,n=p) ).subs(x-x_0,0) + \
Rat(1,2)*( f(x-2*h).series(x-2*h,x0=x_0,n=p)).subs(x-x_0,0) )









    Out[6]:





$$h \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{3} + h^{2} x + h x^{2} + x^{3}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$



In [7]:

    
( Rat(1)*f(x).series(x,x0=x_0,n=p).subs(x-x_0,0) + \
  Rat(-4,3)*( f(x-h).series(x-h,x0=x_0,n=p) ).subs(x-x_0,0) + \
Rat(1,3)*( f(x-2*h).series(x-2*h,x0=x_0,n=p)).subs(x-x_0,0) )









    Out[7]:





$$\frac{2 h}{3} \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(h^{3} + h^{2} x + h x^{2} + x^{3}; \left ( h, \quad x\right )\rightarrow\left ( 0, \quad 0\right )\right)$$

Then "set this" first order partial derivative to zero and then you obtain the "scheme" to "discretize" the boundary condition.

$u_{i,J} = \frac{ 4 u_{i,J-1} - u_{i,J-2} }{ 3 } $

$p_{i,J} = \frac{ 4p_{i,J-1} -p_{i,J-2} }{ 3 } $



In [ ]: