You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

602 lines
20 KiB

5 months ago
from sympy.core import S, Rational
from sympy.integrals.quadrature import (gauss_legendre, gauss_laguerre,
gauss_hermite, gauss_gen_laguerre,
gauss_chebyshev_t, gauss_chebyshev_u,
gauss_jacobi, gauss_lobatto)
def test_legendre():
x, w = gauss_legendre(1, 17)
assert [str(r) for r in x] == ['0']
assert [str(r) for r in w] == ['2.0000000000000000']
x, w = gauss_legendre(2, 17)
assert [str(r) for r in x] == [
'-0.57735026918962576',
'0.57735026918962576']
assert [str(r) for r in w] == [
'1.0000000000000000',
'1.0000000000000000']
x, w = gauss_legendre(3, 17)
assert [str(r) for r in x] == [
'-0.77459666924148338',
'0',
'0.77459666924148338']
assert [str(r) for r in w] == [
'0.55555555555555556',
'0.88888888888888889',
'0.55555555555555556']
x, w = gauss_legendre(4, 17)
assert [str(r) for r in x] == [
'-0.86113631159405258',
'-0.33998104358485626',
'0.33998104358485626',
'0.86113631159405258']
assert [str(r) for r in w] == [
'0.34785484513745386',
'0.65214515486254614',
'0.65214515486254614',
'0.34785484513745386']
def test_legendre_precise():
x, w = gauss_legendre(3, 40)
assert [str(r) for r in x] == [
'-0.7745966692414833770358530799564799221666',
'0',
'0.7745966692414833770358530799564799221666']
assert [str(r) for r in w] == [
'0.5555555555555555555555555555555555555556',
'0.8888888888888888888888888888888888888889',
'0.5555555555555555555555555555555555555556']
def test_laguerre():
x, w = gauss_laguerre(1, 17)
assert [str(r) for r in x] == ['1.0000000000000000']
assert [str(r) for r in w] == ['1.0000000000000000']
x, w = gauss_laguerre(2, 17)
assert [str(r) for r in x] == [
'0.58578643762690495',
'3.4142135623730950']
assert [str(r) for r in w] == [
'0.85355339059327376',
'0.14644660940672624']
x, w = gauss_laguerre(3, 17)
assert [str(r) for r in x] == [
'0.41577455678347908',
'2.2942803602790417',
'6.2899450829374792',
]
assert [str(r) for r in w] == [
'0.71109300992917302',
'0.27851773356924085',
'0.010389256501586136',
]
x, w = gauss_laguerre(4, 17)
assert [str(r) for r in x] == [
'0.32254768961939231',
'1.7457611011583466',
'4.5366202969211280',
'9.3950709123011331']
assert [str(r) for r in w] == [
'0.60315410434163360',
'0.35741869243779969',
'0.038887908515005384',
'0.00053929470556132745']
x, w = gauss_laguerre(5, 17)
assert [str(r) for r in x] == [
'0.26356031971814091',
'1.4134030591065168',
'3.5964257710407221',
'7.0858100058588376',
'12.640800844275783']
assert [str(r) for r in w] == [
'0.52175561058280865',
'0.39866681108317593',
'0.075942449681707595',
'0.0036117586799220485',
'2.3369972385776228e-5']
def test_laguerre_precise():
x, w = gauss_laguerre(3, 40)
assert [str(r) for r in x] == [
'0.4157745567834790833115338731282744735466',
'2.294280360279041719822050361359593868960',
'6.289945082937479196866415765512131657493']
assert [str(r) for r in w] == [
'0.7110930099291730154495901911425944313094',
'0.2785177335692408488014448884567264810349',
'0.01038925650158613574896492040067908765572']
def test_hermite():
x, w = gauss_hermite(1, 17)
assert [str(r) for r in x] == ['0']
assert [str(r) for r in w] == ['1.7724538509055160']
x, w = gauss_hermite(2, 17)
assert [str(r) for r in x] == [
'-0.70710678118654752',
'0.70710678118654752']
assert [str(r) for r in w] == [
'0.88622692545275801',
'0.88622692545275801']
x, w = gauss_hermite(3, 17)
assert [str(r) for r in x] == [
'-1.2247448713915890',
'0',
'1.2247448713915890']
assert [str(r) for r in w] == [
'0.29540897515091934',
'1.1816359006036774',
'0.29540897515091934']
x, w = gauss_hermite(4, 17)
assert [str(r) for r in x] == [
'-1.6506801238857846',
'-0.52464762327529032',
'0.52464762327529032',
'1.6506801238857846']
assert [str(r) for r in w] == [
'0.081312835447245177',
'0.80491409000551284',
'0.80491409000551284',
'0.081312835447245177']
x, w = gauss_hermite(5, 17)
assert [str(r) for r in x] == [
'-2.0201828704560856',
'-0.95857246461381851',
'0',
'0.95857246461381851',
'2.0201828704560856']
assert [str(r) for r in w] == [
'0.019953242059045913',
'0.39361932315224116',
'0.94530872048294188',
'0.39361932315224116',
'0.019953242059045913']
def test_hermite_precise():
x, w = gauss_hermite(3, 40)
assert [str(r) for r in x] == [
'-1.224744871391589049098642037352945695983',
'0',
'1.224744871391589049098642037352945695983']
assert [str(r) for r in w] == [
'0.2954089751509193378830279138901908637996',
'1.181635900603677351532111655560763455198',
'0.2954089751509193378830279138901908637996']
def test_gen_laguerre():
x, w = gauss_gen_laguerre(1, Rational(-1, 2), 17)
assert [str(r) for r in x] == ['0.50000000000000000']
assert [str(r) for r in w] == ['1.7724538509055160']
x, w = gauss_gen_laguerre(2, Rational(-1, 2), 17)
assert [str(r) for r in x] == [
'0.27525512860841095',
'2.7247448713915890']
assert [str(r) for r in w] == [
'1.6098281800110257',
'0.16262567089449035']
x, w = gauss_gen_laguerre(3, Rational(-1, 2), 17)
assert [str(r) for r in x] == [
'0.19016350919348813',
'1.7844927485432516',
'5.5253437422632603']
assert [str(r) for r in w] == [
'1.4492591904487850',
'0.31413464064571329',
'0.0090600198110176913']
x, w = gauss_gen_laguerre(4, Rational(-1, 2), 17)
assert [str(r) for r in x] == [
'0.14530352150331709',
'1.3390972881263614',
'3.9269635013582872',
'8.5886356890120343']
assert [str(r) for r in w] == [
'1.3222940251164826',
'0.41560465162978376',
'0.034155966014826951',
'0.00039920814442273524']
x, w = gauss_gen_laguerre(5, Rational(-1, 2), 17)
assert [str(r) for r in x] == [
'0.11758132021177814',
'1.0745620124369040',
'3.0859374437175500',
'6.4147297336620305',
'11.807189489971737']
assert [str(r) for r in w] == [
'1.2217252674706516',
'0.48027722216462937',
'0.067748788910962126',
'0.0026872914935624654',
'1.5280865710465241e-5']
x, w = gauss_gen_laguerre(1, 2, 17)
assert [str(r) for r in x] == ['3.0000000000000000']
assert [str(r) for r in w] == ['2.0000000000000000']
x, w = gauss_gen_laguerre(2, 2, 17)
assert [str(r) for r in x] == [
'2.0000000000000000',
'6.0000000000000000']
assert [str(r) for r in w] == [
'1.5000000000000000',
'0.50000000000000000']
x, w = gauss_gen_laguerre(3, 2, 17)
assert [str(r) for r in x] == [
'1.5173870806774125',
'4.3115831337195203',
'9.1710297856030672']
assert [str(r) for r in w] == [
'1.0374949614904253',
'0.90575000470306537',
'0.056755033806509347']
x, w = gauss_gen_laguerre(4, 2, 17)
assert [str(r) for r in x] == [
'1.2267632635003021',
'3.4125073586969460',
'6.9026926058516134',
'12.458036771951139']
assert [str(r) for r in w] == [
'0.72552499769865438',
'1.0634242919791946',
'0.20669613102835355',
'0.0043545792937974889']
x, w = gauss_gen_laguerre(5, 2, 17)
assert [str(r) for r in x] == [
'1.0311091440933816',
'2.8372128239538217',
'5.6202942725987079',
'9.6829098376640271',
'15.828473921690062']
assert [str(r) for r in w] == [
'0.52091739683509184',
'1.0667059331592211',
'0.38354972366693113',
'0.028564233532974658',
'0.00026271280578124935']
def test_gen_laguerre_precise():
x, w = gauss_gen_laguerre(3, Rational(-1, 2), 40)
assert [str(r) for r in x] == [
'0.1901635091934881328718554276203028970878',
'1.784492748543251591186722461957367638500',
'5.525343742263260275941422110422329464413']
assert [str(r) for r in w] == [
'1.449259190448785048183829411195134343108',
'0.3141346406457132878326231270167565378246',
'0.009060019811017691281714945129254301865020']
x, w = gauss_gen_laguerre(3, 2, 40)
assert [str(r) for r in x] == [
'1.517387080677412495020323111016672547482',
'4.311583133719520302881184669723530562299',
'9.171029785603067202098492219259796890218']
assert [str(r) for r in w] == [
'1.037494961490425285817554606541269153041',
'0.9057500047030653669269785048806009945254',
'0.05675503380650934725546688857812985243312']
def test_chebyshev_t():
x, w = gauss_chebyshev_t(1, 17)
assert [str(r) for r in x] == ['0']
assert [str(r) for r in w] == ['3.1415926535897932']
x, w = gauss_chebyshev_t(2, 17)
assert [str(r) for r in x] == [
'0.70710678118654752',
'-0.70710678118654752']
assert [str(r) for r in w] == [
'1.5707963267948966',
'1.5707963267948966']
x, w = gauss_chebyshev_t(3, 17)
assert [str(r) for r in x] == [
'0.86602540378443865',
'0',
'-0.86602540378443865']
assert [str(r) for r in w] == [
'1.0471975511965977',
'1.0471975511965977',
'1.0471975511965977']
x, w = gauss_chebyshev_t(4, 17)
assert [str(r) for r in x] == [
'0.92387953251128676',
'0.38268343236508977',
'-0.38268343236508977',
'-0.92387953251128676']
assert [str(r) for r in w] == [
'0.78539816339744831',
'0.78539816339744831',
'0.78539816339744831',
'0.78539816339744831']
x, w = gauss_chebyshev_t(5, 17)
assert [str(r) for r in x] == [
'0.95105651629515357',
'0.58778525229247313',
'0',
'-0.58778525229247313',
'-0.95105651629515357']
assert [str(r) for r in w] == [
'0.62831853071795865',
'0.62831853071795865',
'0.62831853071795865',
'0.62831853071795865',
'0.62831853071795865']
def test_chebyshev_t_precise():
x, w = gauss_chebyshev_t(3, 40)
assert [str(r) for r in x] == [
'0.8660254037844386467637231707529361834714',
'0',
'-0.8660254037844386467637231707529361834714']
assert [str(r) for r in w] == [
'1.047197551196597746154214461093167628066',
'1.047197551196597746154214461093167628066',
'1.047197551196597746154214461093167628066']
def test_chebyshev_u():
x, w = gauss_chebyshev_u(1, 17)
assert [str(r) for r in x] == ['0']
assert [str(r) for r in w] == ['1.5707963267948966']
x, w = gauss_chebyshev_u(2, 17)
assert [str(r) for r in x] == [
'0.50000000000000000',
'-0.50000000000000000']
assert [str(r) for r in w] == [
'0.78539816339744831',
'0.78539816339744831']
x, w = gauss_chebyshev_u(3, 17)
assert [str(r) for r in x] == [
'0.70710678118654752',
'0',
'-0.70710678118654752']
assert [str(r) for r in w] == [
'0.39269908169872415',
'0.78539816339744831',
'0.39269908169872415']
x, w = gauss_chebyshev_u(4, 17)
assert [str(r) for r in x] == [
'0.80901699437494742',
'0.30901699437494742',
'-0.30901699437494742',
'-0.80901699437494742']
assert [str(r) for r in w] == [
'0.21707871342270599',
'0.56831944997474231',
'0.56831944997474231',
'0.21707871342270599']
x, w = gauss_chebyshev_u(5, 17)
assert [str(r) for r in x] == [
'0.86602540378443865',
'0.50000000000000000',
'0',
'-0.50000000000000000',
'-0.86602540378443865']
assert [str(r) for r in w] == [
'0.13089969389957472',
'0.39269908169872415',
'0.52359877559829887',
'0.39269908169872415',
'0.13089969389957472']
def test_chebyshev_u_precise():
x, w = gauss_chebyshev_u(3, 40)
assert [str(r) for r in x] == [
'0.7071067811865475244008443621048490392848',
'0',
'-0.7071067811865475244008443621048490392848']
assert [str(r) for r in w] == [
'0.3926990816987241548078304229099378605246',
'0.7853981633974483096156608458198757210493',
'0.3926990816987241548078304229099378605246']
def test_jacobi():
x, w = gauss_jacobi(1, Rational(-1, 2), S.Half, 17)
assert [str(r) for r in x] == ['0.50000000000000000']
assert [str(r) for r in w] == ['3.1415926535897932']
x, w = gauss_jacobi(2, Rational(-1, 2), S.Half, 17)
assert [str(r) for r in x] == [
'-0.30901699437494742',
'0.80901699437494742']
assert [str(r) for r in w] == [
'0.86831485369082398',
'2.2732777998989693']
x, w = gauss_jacobi(3, Rational(-1, 2), S.Half, 17)
assert [str(r) for r in x] == [
'-0.62348980185873353',
'0.22252093395631440',
'0.90096886790241913']
assert [str(r) for r in w] == [
'0.33795476356635433',
'1.0973322242791115',
'1.7063056657443274']
x, w = gauss_jacobi(4, Rational(-1, 2), S.Half, 17)
assert [str(r) for r in x] == [
'-0.76604444311897804',
'-0.17364817766693035',
'0.50000000000000000',
'0.93969262078590838']
assert [str(r) for r in w] == [
'0.16333179083642836',
'0.57690240318269103',
'1.0471975511965977',
'1.3541609083740761']
x, w = gauss_jacobi(5, Rational(-1, 2), S.Half, 17)
assert [str(r) for r in x] == [
'-0.84125353283118117',
'-0.41541501300188643',
'0.14231483827328514',
'0.65486073394528506',
'0.95949297361449739']
assert [str(r) for r in w] == [
'0.090675770007435372',
'0.33391416373675607',
'0.65248870981926643',
'0.94525424081394926',
'1.1192597692123861']
x, w = gauss_jacobi(1, 2, 3, 17)
assert [str(r) for r in x] == ['0.14285714285714286']
assert [str(r) for r in w] == ['1.0666666666666667']
x, w = gauss_jacobi(2, 2, 3, 17)
assert [str(r) for r in x] == [
'-0.24025307335204215',
'0.46247529557426437']
assert [str(r) for r in w] == [
'0.48514624517838660',
'0.58152042148828007']
x, w = gauss_jacobi(3, 2, 3, 17)
assert [str(r) for r in x] == [
'-0.46115870378089762',
'0.10438533038323902',
'0.62950064612493132']
assert [str(r) for r in w] == [
'0.17937613502213266',
'0.61595640991147154',
'0.27133412173306246']
x, w = gauss_jacobi(4, 2, 3, 17)
assert [str(r) for r in x] == [
'-0.59903470850824782',
'-0.14761105199952565',
'0.32554377081188859',
'0.72879429738819258']
assert [str(r) for r in w] == [
'0.067809641836772187',
'0.38956404952032481',
'0.47995970868024150',
'0.12933326662932816']
x, w = gauss_jacobi(5, 2, 3, 17)
assert [str(r) for r in x] == [
'-0.69045775012676106',
'-0.32651993134900065',
'0.082337849552034905',
'0.47517887061283164',
'0.79279429464422850']
assert [str(r) for r in w] == [
'0.027410178066337099',
'0.21291786060364828',
'0.43908437944395081',
'0.32220656547221822',
'0.065047683080512268']
def test_jacobi_precise():
x, w = gauss_jacobi(3, Rational(-1, 2), S.Half, 40)
assert [str(r) for r in x] == [
'-0.6234898018587335305250048840042398106323',
'0.2225209339563144042889025644967947594664',
'0.9009688679024191262361023195074450511659']
assert [str(r) for r in w] == [
'0.3379547635663543330553835737094171534907',
'1.097332224279111467485302294320899710461',
'1.706305665744327437921957515249186020246']
x, w = gauss_jacobi(3, 2, 3, 40)
assert [str(r) for r in x] == [
'-0.4611587037808976179121958105554375981274',
'0.1043853303832390210914918407615869143233',
'0.6295006461249313240934312425211234110769']
assert [str(r) for r in w] == [
'0.1793761350221326596137764371503859752628',
'0.6159564099114715430909548532229749439714',
'0.2713341217330624639619353762933057474325']
def test_lobatto():
x, w = gauss_lobatto(2, 17)
assert [str(r) for r in x] == [
'-1',
'1']
assert [str(r) for r in w] == [
'1.0000000000000000',
'1.0000000000000000']
x, w = gauss_lobatto(3, 17)
assert [str(r) for r in x] == [
'-1',
'0',
'1']
assert [str(r) for r in w] == [
'0.33333333333333333',
'1.3333333333333333',
'0.33333333333333333']
x, w = gauss_lobatto(4, 17)
assert [str(r) for r in x] == [
'-1',
'-0.44721359549995794',
'0.44721359549995794',
'1']
assert [str(r) for r in w] == [
'0.16666666666666667',
'0.83333333333333333',
'0.83333333333333333',
'0.16666666666666667']
x, w = gauss_lobatto(5, 17)
assert [str(r) for r in x] == [
'-1',
'-0.65465367070797714',
'0',
'0.65465367070797714',
'1']
assert [str(r) for r in w] == [
'0.10000000000000000',
'0.54444444444444444',
'0.71111111111111111',
'0.54444444444444444',
'0.10000000000000000']
def test_lobatto_precise():
x, w = gauss_lobatto(3, 40)
assert [str(r) for r in x] == [
'-1',
'0',
'1']
assert [str(r) for r in w] == [
'0.3333333333333333333333333333333333333333',
'1.333333333333333333333333333333333333333',
'0.3333333333333333333333333333333333333333']