CoCalc -- volume_integral

| Hosted by CoCalc | Download

import numpy as np
%typeset_mode False

var('x, xi')

x_nodes = [-1, 1]
xi_LGL = [-1.0,           \
          -0.87174014851, \
          -0.591700181433,\
          -0.209299217902,\
          0.209299217902, \
          0.591700181433, \
          0.87174014851,  \
          1.0]

lobatto_weights = [0.03571428571429,\
                   0.21070422714351,\
                   0.34112269248350,\
                   0.41245879465870,\
                   0.41245879465870,\
                   0.34112269248350,\
                   0.21070422714351,\
                   0.03571428571429]

(x, xi)

numerical_integral(np.e ** (- x**2/(0.4 ** 2)), -1, 1)

(0.7086930188940249, 7.868073069672576e-15)

    #Assigning L_0 and L_1 as the lagrange basis polynomials.
L = [0]
L_0 =     ((xi - (xi_LGL[1])) / (xi_LGL[0] - xi_LGL[1]))\
        * ((xi - (xi_LGL[2])) / (xi_LGL[0] - xi_LGL[2]))\
        * ((xi - (xi_LGL[3])) / (xi_LGL[0] - xi_LGL[3]))\
        * ((xi - (xi_LGL[4])) / (xi_LGL[0] - xi_LGL[4]))\
        * ((xi - (xi_LGL[5])) / (xi_LGL[0] - xi_LGL[5]))\
        * ((xi - (xi_LGL[6])) / (xi_LGL[0] - xi_LGL[6]))\
        * ((xi - (xi_LGL[7])) / (xi_LGL[0] - xi_LGL[7]))

L_1 =     ((xi - (xi_LGL[0])) / (xi_LGL[1] - xi_LGL[0]))\
        * ((xi - (xi_LGL[2])) / (xi_LGL[1] - xi_LGL[2]))\
        * ((xi - (xi_LGL[3])) / (xi_LGL[1] - xi_LGL[3]))\
        * ((xi - (xi_LGL[4])) / (xi_LGL[1] - xi_LGL[4]))\
        * ((xi - (xi_LGL[5])) / (xi_LGL[1] - xi_LGL[5]))\
        * ((xi - (xi_LGL[6])) / (xi_LGL[1] - xi_LGL[6]))\
        * ((xi - (xi_LGL[7])) / (xi_LGL[1] - xi_LGL[7]))
        
L_2 =   ((xi - (xi_LGL[0])) / (xi_LGL[2] - xi_LGL[0])) \
      * ((xi - (xi_LGL[1])) / (xi_LGL[2] - xi_LGL[1])) \
      * ((xi - (xi_LGL[3])) / (xi_LGL[2] - xi_LGL[3])) \
      * ((xi - (xi_LGL[4])) / (xi_LGL[2] - xi_LGL[4])) \
      * ((xi - (xi_LGL[5])) / (xi_LGL[2] - xi_LGL[5])) \
      * ((xi - (xi_LGL[6])) / (xi_LGL[2] - xi_LGL[6])) \
      * ((xi - (xi_LGL[7])) / (xi_LGL[2] - xi_LGL[7]))

L_3 =   ((xi - (xi_LGL[0])) / (xi_LGL[3] - xi_LGL[0])) \
      * ((xi - (xi_LGL[1])) / (xi_LGL[3] - xi_LGL[1])) \
      * ((xi - (xi_LGL[2])) / (xi_LGL[3] - xi_LGL[2])) \
      * ((xi - (xi_LGL[4])) / (xi_LGL[3] - xi_LGL[4])) \
      * ((xi - (xi_LGL[5])) / (xi_LGL[3] - xi_LGL[5])) \
      * ((xi - (xi_LGL[6])) / (xi_LGL[3] - xi_LGL[6])) \
      * ((xi - (xi_LGL[7])) / (xi_LGL[3] - xi_LGL[7]))
        
L_4 =   ((xi - (xi_LGL[0])) / (xi_LGL[4] - xi_LGL[0])) \
      * ((xi - (xi_LGL[1])) / (xi_LGL[4] - xi_LGL[1])) \
      * ((xi - (xi_LGL[2])) / (xi_LGL[4] - xi_LGL[2])) \
      * ((xi - (xi_LGL[3])) / (xi_LGL[4] - xi_LGL[3])) \
      * ((xi - (xi_LGL[5])) / (xi_LGL[4] - xi_LGL[5])) \
      * ((xi - (xi_LGL[6])) / (xi_LGL[4] - xi_LGL[6])) \
      * ((xi - (xi_LGL[7])) / (xi_LGL[4] - xi_LGL[7]))


L_5 =   ((xi - (xi_LGL[0])) / (xi_LGL[5] - xi_LGL[0])) \
      * ((xi - (xi_LGL[1])) / (xi_LGL[5] - xi_LGL[1])) \
      * ((xi - (xi_LGL[2])) / (xi_LGL[5] - xi_LGL[2])) \
      * ((xi - (xi_LGL[3])) / (xi_LGL[5] - xi_LGL[3])) \
      * ((xi - (xi_LGL[4])) / (xi_LGL[5] - xi_LGL[4])) \
      * ((xi - (xi_LGL[6])) / (xi_LGL[5] - xi_LGL[6])) \
      * ((xi - (xi_LGL[7])) / (xi_LGL[5] - xi_LGL[7]))
        
L_6 =   ((xi - (xi_LGL[0])) / (xi_LGL[6] - xi_LGL[0])) \
      * ((xi - (xi_LGL[1])) / (xi_LGL[6] - xi_LGL[1])) \
      * ((xi - (xi_LGL[2])) / (xi_LGL[6] - xi_LGL[2])) \
      * ((xi - (xi_LGL[3])) / (xi_LGL[6] - xi_LGL[3])) \
      * ((xi - (xi_LGL[4])) / (xi_LGL[6] - xi_LGL[4])) \
      * ((xi - (xi_LGL[5])) / (xi_LGL[6] - xi_LGL[5])) \
      * ((xi - (xi_LGL[7])) / (xi_LGL[6] - xi_LGL[7]))
        

L_7 =   ((xi - (xi_LGL[0])) / (xi_LGL[7] - xi_LGL[0])) \
      * ((xi - (xi_LGL[1])) / (xi_LGL[7] - xi_LGL[1])) \
      * ((xi - (xi_LGL[2])) / (xi_LGL[7] - xi_LGL[2])) \
      * ((xi - (xi_LGL[3])) / (xi_LGL[7] - xi_LGL[3])) \
      * ((xi - (xi_LGL[4])) / (xi_LGL[7] - xi_LGL[4])) \
      * ((xi - (xi_LGL[5])) / (xi_LGL[7] - xi_LGL[5])) \
      * ((xi - (xi_LGL[6])) / (xi_LGL[7] - xi_LGL[6])) \

print(L[0]    )

0

print(e ** (- (xi_LGL * xi_LGL) / (0.4 ** 2)))

Error in lines 1-1
Traceback (most recent call last):
  File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 996, in execute
    exec compile(block+'\n', '', 'single') in namespace, locals
  File "", line 1, in <module>
TypeError: can't multiply sequence by non-int of type 'list'

flux = e ** (- (xi ** 2) / (0.4 ** 2))
print(numerical_integral(diff(L_0) * flux, -1, 1)[0])
print(numerical_integral(diff(L_1) * flux, -1, 1)[0])
print(numerical_integral(diff(L_2) * flux, -1, 1)[0])
print(numerical_integral(diff(L_3) * flux, -1, 1)[0])
print(numerical_integral(diff(L_4) * flux, -1, 1)[0])
print(numerical_integral(diff(L_5) * flux, -1, 1)[0])
print(numerical_integral(diff(L_6) * flux, -1, 1)[0])
print(numerical_integral(diff(L_7) * flux, -1, 1)[0])

-0.0118823589586
0.0086260537921
-0.3480503904
-0.688212551383
0.688212551383
0.3480503904
-0.0086260537921
0.0118823589586


for i in range(10):
    element_limits = [-1 + 0.2 * i, -0.8 + 0.2 * i]
    x_element = sum(element_limits) / 2 + xi * (element_limits[1] - element_limits[0]) / 2
    flux = e ** (- (x_element ** 2) / (0.4 ** 2))
    print(numerical_integral(diff(L_0) * flux, -1, 1)[0])
    print(numerical_integral(diff(L_1) * flux, -1, 1)[0])
    print(numerical_integral(diff(L_2) * flux, -1, 1)[0])
    print(numerical_integral(diff(L_3) * flux, -1, 1)[0])
    print(numerical_integral(diff(L_4) * flux, -1, 1)[0])
    print(numerical_integral(diff(L_5) * flux, -1, 1)[0])
    print(numerical_integral(diff(L_6) * flux, -1, 1)[0])
    print(numerical_integral(diff(L_7) * flux, -1, 1)[0])
    print('--------------')

-0.00201663487667
-0.000588597708116
-0.00130167737191
-0.00236838757932
-0.00362050204766
-0.00432019709409
-0.00344551201015
0.0176615086879
--------------
-0.018969769374
-0.00431252844519
-0.00882630935977
-0.0144355176966
-0.019612124119
-0.0209837936827
-0.0154359890788
0.102576031756
--------------
-0.108222418798
-0.0179274222595
-0.0337807018822
-0.0492589052599
-0.0588472807471
-0.0557970236273
-0.0374764132459
0.361310165819
--------------
-0.374448714304
-0.0399576371245
-0.0683852285846
-0.0869229749357
-0.0884322503841
-0.0714664112839
-0.0422339853622
0.771847201979
--------------
-0.785754362849
-0.0396035640187
-0.0579313769517
-0.0569022801117
-0.0392041960688
-0.0172295769141
-0.00337464521455
1.00000000213
--------------
-1.00000000213
0.00337464521455
0.0172295769141
0.0392041960688
0.0569022801117
0.0579313769517
0.0396035640187
0.785754362849
--------------
-0.771847201979
0.0422339853622
0.0714664112839
0.0884322503841
0.0869229749357
0.0683852285846
0.0399576371245
0.374448714304
--------------
-0.361310165819
0.0374764132459
0.0557970236273
0.0588472807471
0.0492589052599
0.0337807018822
0.0179274222595
0.108222418798
--------------
-0.102576031756
0.0154359890788
0.0209837936827
0.019612124119
0.0144355176966
0.00882630935977
0.00431252844519
0.018969769374
--------------
-0.0176615086879
0.00344551201015
0.00432019709409
0.00362050204766
0.00236838757932
0.00130167737191
0.000588597708116
0.00201663487667
--------------

-0.00201663487667 - -0.00201663465967866

-2.16991339827000e-10

#Calculating the flux integral using Gauss-Lobatto quadrature
lobatto_integral = lobatto_weights[0] * diff(L_0, xi)(xi_LGL[0]) * flux(element_0_nodes[0])\
                  +lobatto_weights[1] * diff(L_0, xi)(xi_LGL[1]) * flux(element_0_nodes[1])\
                  +lobatto_weights[2] * diff(L_0, xi)(xi_LGL[2]) * flux(element_0_nodes[2])\
                  +lobatto_weights[3] * diff(L_0, xi)(xi_LGL[3]) * flux(element_0_nodes[3])\
                  +lobatto_weights[4] * diff(L_0, xi)(xi_LGL[4]) * flux(element_0_nodes[4])\
                  +lobatto_weights[5] * diff(L_0, xi)(xi_LGL[5]) * flux(element_0_nodes[5])\
                  +lobatto_weights[6] * diff(L_0, xi)(xi_LGL[6]) * flux(element_0_nodes[6])\
                  +lobatto_weights[7] * diff(L_0, xi)(xi_LGL[7]) * flux(element_0_nodes[7])

print(lobatto_integral)

-0.00201663465967866

f = xi ** 2 + xi ** 3

f(2 * xi)

\displaystyle 8 \, \xi^{3} + 4 \, \xi^{2}


diff(L_0, xi)(xi_LGL[0])
diff(L_0, xi)(xi_LGL[1])
diff(L_0, xi)(xi_LGL[2])
diff(L_0, xi)(xi_LGL[3])
diff(L_0, xi)(xi_LGL[4])
diff(L_0, xi)(xi_LGL[5])
diff(L_0, xi)(xi_LGL[6])
diff(L_0, xi)(xi_LGL[7])

-14.0000000000226
-3.20991570302344
0.792476681323880
-0.372150435728984
0.243330712724289
-0.203284568901545
0.219957514771985
-0.500000000000000