CoCalc Public FilesPM_2_5 / wave_equation / worksheets / volume_integral_flux.sagewsOpen in with one click!
Authors: Mani Chandra, Balavarun P, Aman Abhishek Tiwari
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)
8ξ3+4ξ2\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