CoCalc Public FilesPM_2_5 / wave_equation / worksheets / unitTestA_Matrix.sagews
%typeset_mode False
from scipy import special as sp
from scipy import integrate
from sage.symbolic.integration.integral import indefinite_integral
import numpy as np
import math


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)
L = [xi, xi, xi, xi, xi, xi, xi, xi]

#Assigning L_0 and L_1 as the lagrange basis polynomials.
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]))




numerical_integral(L[0] * L[1], -1, 1)[0]

0.005783175201965206



for i in range (0, 8):
print(i, '\n')
for j in range(0, 8):
numerical_integral(L[i] * L[j], -1, 1)[0]

(0, '\n') 0.03333333333332194 0.005783175201965206 -0.007358427761753982 0.008091331778355441 -0.008091331778233877 0.007358427761705623 -0.005783175202249493 0.002380952380963754 (1, '\n') 0.005783175201965206 0.19665727866729804 0.017873132323192046 -0.01965330750343234 0.019653307503020866 -0.017873132322725152 0.014046948476303067 -0.005783175202249493 (2, '\n') -0.007358427761753982 0.017873132323192046 0.31838117965137114 0.025006581762566073 -0.025006581761945083 0.022741512832051156 -0.017873132322725152 0.007358427761705624 (3, '\n') 0.008091331778355441 -0.01965330750343234 0.025006581762566073 0.3849615416814164 0.027497252976343693 -0.025006581761945083 0.019653307503020863 -0.008091331778233875 (4, '\n') -0.008091331778233877 0.019653307503020866 -0.025006581761945083 0.027497252976343693 0.3849615416814164 0.025006581762566073 -0.019653307503432346 0.008091331778355443 (5, '\n') 0.007358427761705623 -0.017873132322725152 0.022741512832051156 -0.025006581761945083 0.025006581762566073 0.31838117965137114 0.017873132323192046 -0.0073584277617539835 (6, '\n') -0.005783175202249493 0.014046948476303067 -0.017873132322725152 0.019653307503020863 -0.019653307503432346 0.017873132323192046 0.19665727866729804 0.0057831752019652065 (7, '\n') 0.002380952380963754 -0.005783175202249493 0.007358427761705624 -0.008091331778233875 0.008091331778355443 -0.0073584277617539835 0.0057831752019652065 0.033333333333321946


LGL_list = [ \
[-1.0,1.0],                                                               \
[-1.0,0.0,1.0],                                                           \
[-1.0,-0.4472135955,0.4472135955,1.0],                                    \
[-1.0,-0.654653670708,0.0,0.654653670708,1.0],                            \
[-1.0,-0.765055323929,-0.285231516481,0.285231516481,0.765055323929,1.0], \
[-1.0,-0.830223896279,-0.468848793471,0.0,0.468848793471,0.830223896279,  \
1.0],                                                                     \
[-1.0,-0.87174014851,-0.591700181433,-0.209299217902,0.209299217902,      \
0.591700181433,0.87174014851,1.0],                                        \
[-1.0,-0.899757995411,-0.677186279511,-0.363117463826,0.0,0.363117463826, \
0.677186279511,0.899757995411,1.0],                                       \
[-1.0,-0.919533908167,-0.738773865105,-0.47792494981,-0.165278957666,     \
0.165278957666,0.47792494981,0.738773865106,0.919533908166,1.0],          \
[-1.0,-0.934001430408,-0.784483473663,-0.565235326996,-0.295758135587,    \
0.0,0.295758135587,0.565235326996,0.784483473663,0.934001430408,1.0],     \
[-1.0,-0.944899272223,-0.819279321644,-0.632876153032,-0.399530940965,    \
-0.136552932855,0.136552932855,0.399530940965,0.632876153032,             \
0.819279321644,0.944899272223,1.0],                                       \
[-1.0,-0.953309846642,-0.846347564652,-0.686188469082,-0.482909821091,    \
-0.249286930106,0.0,0.249286930106,0.482909821091,0.686188469082,         \
0.846347564652,0.953309846642,1.0],                                       \
[-0.999999999996,-0.959935045274,-0.867801053826,-0.728868599093,         \
-0.550639402928,-0.342724013343,-0.116331868884,0.116331868884,           \
0.342724013343,0.550639402929,0.728868599091,0.86780105383,               \
0.959935045267,1.0],                                                      \
[-0.999999999996,-0.965245926511,-0.885082044219,-0.763519689953,         \
-0.60625320547,-0.420638054714,-0.215353955364,0.0,0.215353955364,        \
0.420638054714,0.60625320547,0.763519689952,0.885082044223,               \
0.965245926503,1.0],                                                      \
[-0.999999999984,-0.9695680463,-0.899200533072,-0.792008291871,           \
-0.65238870288,-0.486059421887,-0.299830468901,-0.101326273522,           \
0.101326273522,0.299830468901,0.486059421887,0.652388702882,              \
0.792008291863,0.899200533092,0.969568046272,0.999999999999]]             \


def L0_L1(xi_LGL):
xi_LGL = np.array(xi_LGL)
lagrange_basis_poly = []

for j in np.arange(xi_LGL.shape[0]):
lagrange_basis_j = np.poly1d([1])

for m in np.arange(xi_LGL.shape[0]):
if m != j:
lagrange_basis_j *= np.poly1d([1, -xi_LGL[m]]) \
/ (xi_LGL[j] - xi_LGL[m])
lagrange_basis_poly.append(lagrange_basis_j)

return lagrange_basis_poly[0] * lagrange_basis_poly[1]

print(L0_L1(LGL_points))

Error in lines 1-1 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 995, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> File "", line 4, in L0_L1 IndexError: tuple index out of range
numerical_integral(L0_L1(LGL_list[6]), -1, 1)

(0.005783175201964073, 4.71163122450593e-16)
for i in range(0,28):
print(numerical_integral(L0_L1(LGL_list[i]), -1, 1)[0])

0.333333333333 0.133333333333 0.05323971375 0.0259259259259 0.014440434437 0.00883182542957 0.00578317520196 0.00398834520818 0.00286460651796 0.00212583206228 0.00162049687462 0.00126330203786 0.00100374856554 0.000810649719457 0.000664040592284
Error in lines 1-2 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 995, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 2, in <module> IndexError: list index out of range
def LGL_points(N):
xi                 = np.poly1d([1, 0])
legendre_N_minus_1 = N * (xi * sp.legendre(N - 1) - sp.legendre(N))
lgl_points         = legendre_N_minus_1.r
lgl_points.sort()
lgl_points         = (lgl_points)

return lgl_points

print(LGL_points(7))

[-1. -0.8302239 -0.46884879 0. 0.46884879 0.8302239 1. ]
print(numerical_integral(L0_L1(LGL_points(8)), -1, 1))

(0.00578317520220532, 4.711631224498659e-16)
for i in range (3, 30):
print(numerical_integral(L0_L1(LGL_points(i)), -1, 1))

(0.13333333333333341, 2.7731031508685677e-15) (0.05323971374999505, 1.7332838543151235e-15) (0.025925925925926265, 1.1568373972098388e-15) (0.014440434436697796, 8.180747765513089e-16) (0.008831825429837327, 6.058503004859771e-16) (0.00578317520220532, 4.711631224498659e-16) (0.003988345207907132, 3.7296603402372665e-16) (0.0028646065180825096, 3.0374974956255124e-16) (0.0021258320621168345, 2.5226541016198267e-16) (0.0016204968750730782, 5.567157207471083e-16) (0.0012633020355843898, 1.8164735277876645e-16) (0.001003748564545736, 2.2338228618937467e-14) (0.0008106497277544586, 6.454994913541323e-13) (0.0006640406604325044, 1.5566421968657747e-12) (0.0005507393097814919, 4.7083985358726706e-11) (0.0004618008552699972, 2.1039684070047994e-10) (0.00039103100441685167, 7.277576934275988e-10) (0.00033398425097875384, 4.474912702644648e-09) (0.00028756185227232593, 5.894733002978138e-07) (0.00024821818705448266, 6.75715008280175e-07) (0.00021808168014934384, 1.5054891962539568e-06) (0.00019254567973207468, 3.4873510873850225e-05) (0.00019143453010449118, 0.00019140574954508706) (-0.000672459855719221, 0.000349735858050207) (0.001327500977237568, 0.00195693446895412) (0.008536201682912633, 0.009403991970390362) (0.032276067956334815, 0.03684427040090083)
print((L0_L1(LGL_points(30))))

58 57 56 55 -9.096e+12 x + 1.812e+13 x + 1.118e+14 x - 2.406e+14 x 54 53 52 51 - 6.386e+14 x + 1.51e+15 x + 2.241e+15 x - 5.962e+15 x 50 49 48 47 - 5.366e+15 x + 1.66e+16 x + 9.147e+15 x - 3.468e+16 x 46 45 44 43 - 1.107e+16 x + 5.644e+16 x + 8.721e+15 x - 7.334e+16 x 42 41 40 39 - 2.335e+15 x + 7.738e+16 x - 4.86e+15 x - 6.706e+16 x 38 37 36 35 + 9.24e+15 x + 4.81e+16 x - 9.544e+15 x - 2.869e+16 x 34 33 32 31 + 7.122e+15 x + 1.427e+16 x - 4.133e+15 x - 5.914e+15 x 30 29 28 27 + 1.92e+15 x + 2.041e+15 x - 7.23e+14 x - 5.843e+14 x 26 25 24 23 + 2.217e+14 x + 1.38e+14 x - 5.531e+13 x - 2.665e+13 x 22 21 20 19 + 1.118e+13 x + 4.169e+12 x - 1.814e+12 x - 5.209e+11 x 18 17 16 15 + 2.337e+11 x + 5.106e+10 x - 2.35e+10 x - 3.837e+09 x 14 13 12 11 + 1.803e+09 x + 2.143e+08 x - 1.024e+08 x - 8.537e+06 x 10 9 8 7 6 + 4.136e+06 x + 2.285e+05 x - 1.119e+05 x - 3761 x + 1856 x 5 4 3 2 + 32.85 x - 16.3 x - 0.1139 x + 0.05663 x + 0.0001312 x - 6.532e-05
numerical_integral(9.096e12 * x**58, -1, 1)

(308338983050.84705, 0.0034232503837255627)
from sage.symbolic.integration.integral import definite_integral

definite_integral(9.096e12 * x**58,x,-1,1)

308338983050.8475
x = np.poly1d([1, 0])
y = np.poly1d([1, 0])

for i in range(0,57):
x = x*y

print(x)

58 1 x
print(np.poly1d.integ(9.096e12 * x)(1) - np.poly1d.integ(9.096e12 * x)(-1) - 308338983050.)

0.847473144531
definite_integral(9.096e12 * x**58,x,-1,1) - numerical_integral(9.096e12 * x**58, -1, 1)[0]

0.00042724609375
for i in range (3, 30):
print(definite_integral(L0_L1(LGL_points(i)), x, -1, 1))

Error in lines 1-2 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 995, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 2, in <module> File "sage/symbolic/function.pyx", line 996, in sage.symbolic.function.BuiltinFunction.__call__ (/ext/sage/sage-8.0/src/build/cythonized/sage/symbolic/function.cpp:11397) res = super(BuiltinFunction, self).__call__( File "sage/symbolic/function.pyx", line 474, in sage.symbolic.function.Function.__call__ (/ext/sage/sage-8.0/src/build/cythonized/sage/symbolic/function.cpp:6282) raise TypeError("cannot coerce arguments: %s" % (err)) TypeError: cannot coerce arguments: no canonical coercion from <class 'numpy.lib.polynomial.poly1d'> to Symbolic Ring
print(np.poly1d.integ(L0_L1(LGL_points(5))))

9 8 7 6 5 4 -0.1985 x + 0.3695 x + 0.1975 x - 0.7038 x + 0.181 x + 0.3167 x 3 - 0.1671 x
print(L0_L1(LGL_points(5)))

8 7 6 5 4 3 2 -1.786 x + 2.956 x + 1.383 x - 4.223 x + 0.9051 x + 1.267 x - 0.5012 x
for i in range (3, 30):
print(np.poly1d.integ(L0_L1(LGL_points(i)))(1) - np.poly1d.integ(L0_L1(LGL_points(i)))(-1))

0.133333333333 0.05323971375 0.0259259259259 0.0144404344367 0.00883182542984 0.0057831752022 0.00398834520791 0.00286460651807 0.00212583206209 0.00162049687496 0.00126330203518 0.00100374856552 0.000810649721334 0.000664040648327 0.000550739482146 0.000461799455741 0.000391029750504 0.000333980074166 0.000287356982044 0.000250148944115 0.000218478781885 0.00019966664253 0.000120453264427 -0.000364078311827 0.00218563723183 -0.00171491554656 -0.00311789820183