InputMode:=Word

InputMode:=Word

PrecisionDigits:=64

PrecisionDigits:=64

LOAD("c1999.mth")

true

"This is an example of a elliptic integral."

int1:=INT(1/(SQRT((t+a1)*(t+a2)*(t+a3)*(t+a4))*(t+a5)^2*(t+a6)^2),t,0,1)

"This elliptic integral is to this type:"

IA(-2*e5-2*e6)

IA(-2*e5-2*e6)

at:=[a1,a2,a3,a4,a5,a6]

at:=[a1,a2,a3,a4,a5,a6]

bt:=[1,1,1,1,1,1]

bt:=[1,1,1,1,1,1]

mm:=[0,0,0,0,-2,-2]

mm:=[0,0,0,0,-2,-2]

"Since all the values of b are 1, we can simplify the dij difference."

DFA(i,j):=at SUB i-at SUB j

DFA(i,j):=[a1,a2,a3,a4,a5,a6] SUB i-[a1,a2,a3,a4,a5,a6] SUB j

"This is the symbolic solution in terms of basic integrals."

int3:=D219A(1,mm,6,4)

int3:=(2*(a1-a6)*(a2-a6)*(a3-a6)*(a4-a6)*(a5-a6)*AA(e1+e2+e3+e4-e5)+2*(a1-a5)~
*(a2-a5)*(a3-a5)*(a4-a5)*(a5-a6)*AA(e1+e2+e3+e4-e6)+2*(a6-a5)*(a1*(a2*(a3*(2*~
a4-a5-a6)-a4*(a5+a6)+a5^2+a6^2)-a3*(a4*(a5+a6)-a5^2-a6^2)+a4*(a5^2+a6^2)-a5^3~
-a6^3)-a2*(a3*(a4*(a5+a6)-a5^2-a6^2)-a4*(a5^2+a6^2)+a5^3+a6^3)+a3*(a4*(a5^2+a~
6^2)-a5^3-a6^3)-a4*(a5^3+a6^3)+a5^4+a6^4)*AA(e2+e3+e4)+(a1-a2)*(a1-a3)*(a1-a4~
)*(a5-a6)*(a1*(a2*(a3*(2*a4-a5-a6)-a4*(a5+a6)+a5^2+a6^2)-a3*(a4*(a5+a6)-a5^2-~
a6^2)+a4*(a5^2+a6^2)-a5^3-a6^3)-a2*(a3*(a4*(a5+a6)-a5^2-a6^2)-a4*(a5^2+a6^2)+~
a5^3+a6^3)+a3*(a4*(a5^2+a6^2)-a5^3-a6^3)-a4*(a5^3+a6^3)+a5^4+a6^4)*IA(-e1)+(a~
1-a6)*(a2-a6)*(a3-a6)*(a6-a4)*(a1*(a2*(a3*(4*a4-5*a5+a6)+a4*(a6-5*a5)+2*a5*(3~
*a5-a6))+a3*(2*a5*(3*a5-a6)-a4*(5*a5-a6))+a5*(2*a4*(3*a5-a6)-a5*(7*a5-3*a6)))~
+a2*(a5*(2*a4*(3*a5-a6)-a5*(7*a5-3*a6))-a3*(a4*(5*a5-a6)-2*a5*(3*a5-a6)))+a5*~
(a3*(2*a4*(3*a5-a6)-a5*(7*a5-3*a6))-a5*(a4*(7*a5-3*a6)-4*a5*(2*a5-a6))))*IA(-~
e5)+(a1-a5)*((a2-a5)*(a3-a5)*(a4-a5)*(a1*(a2*(a3*(4*a4+a5-5*a6)+a4*(a5-5*a6)-~
2*a6*(a5-3*a6))+a3*(a4*(a5-5*a6)-2*a6*(a5-3*a6))-a6*(2*a4*(a5-3*a6)-a6*(3*a5-~
7*a6)))+a2*(a3*(a4*(a5-5*a6)-2*a6*(a5-3*a6))-a6*(2*a4*(a5-3*a6)-a6*(3*a5-7*a6~
)))-a6*(a3*(2*a4*(a5-3*a6)-a6*(3*a5-7*a6))-a6*(a4*(3*a5-7*a6)-4*a6*(a5-2*a6))~
))*IA(-e6)+(a1-a6)*(a6-a5)*(a1*(a2*(a3*(2*a4-a5-a6)-a4*(a5+a6)+a5^2+a6^2)-a3*~
(a4*(a5+a6)-a5^2-a6^2)+a4*(a5^2+a6^2)-a5^3-a6^3)-a2^2*(a3*(2*a4-a5-a6)-a4*(a5~
+a6)+a5^2+a6^2)-a2*(a3^2*(2*a4-a5-a6)+a3*(a4-2*(a5+a6))*(2*a4-a5-a6)-(a4^2-2*~
a4*(a5+a6)+a5^2+a5*a6+a6^2)*(a5+a6))+a3^2*(a4*(a5+a6)-a5^2-a6^2)+a3*(a4^2-2*a~
4*(a5+a6)+a5^2+a5*a6+a6^2)*(a5+a6)-a4^2*(a5^2+a6^2)+a4*(a5^3+2*a5^2*a6+2*a5*a~
6^2+a6^3)-2*a5*a6*(a5^2+a6^2))*IA(0)))/(2*(a1-a5)*(a1-a6)*(a2-a5)*(a2-a6)*(a3~
-a5)*(a3-a6)*(a4-a5)*(a4-a6)*(a6-a5)^3)

"Check the answer using numerical integration and a random example."

"Generate random values of a."

an:=VECTOR(RANDOM(128)/RANDOM(128),i,1,6)

an:=[43/75,29/65,119/124,73/117,116/69,53/28]

"Substute the values of a in the integral and evaluate it to 12 digits."

int2:=APPROX(SUBST(int1,[a1,a2,a3,a4,a5,a6],an),12)

int2:=1293087/23740243

"Substute the values of a in the symbolic solution."

int4:=SUBST(int3,[a1,a2,a3,a4,a5,a6],an)

int4:=-4232*(1789829891020668558197892015000*AA(e1+e2+e3+e4-e5)+8263198943538~
04850723811360000*AA(e1+e2+e3+e4-e6)-2616149785374473408921703375000*AA(e2+e3~
+e4)+3252069309886469960574103408*IA(-e1)+15*(924630970501288612243636805979*~
IA(-e5)-637*(2170247831638616325279224758*IA(-e6)-330652465111719637989415105~
*IA(0))))/354261834356619312270486986666175

"Expand the symbolic solution to reveal the basic integrals."

int5:=TERMS(EXPAND(int4))

int5:=[-69585631554628800*AA(e1+e2+e3+e4-e5)/3254516850861583,-50082214766438~
400*AA(e1+e2+e3+e4-e6)/5073574307278031,3084876849714687749721480000*AA(e2+e3~
+e4)/98708359449876165832998833,-33649773397162691621392678784*IA(-e1)/866165~
854172663355184564759575,-14335213855710841650792*IA(-e5)/8652133048015518405~
5,23131892295883231600208*IA(-e6)/93379135125452160555,-120999525094057608342~
40952*IA(0)/320597342132401837783107]

"The first term is not a basic integral, but a function that can be evaluated~
."

int5 SUB 1

-69585631554628800*AA(e1+e2+e3+e4-e5)/3254516850861583

"The function is:"

AA(e1+e2+e3+e4-e5)=AF([1,1,1,1,-1,0],4,6,an,bt,1,0)

AA(e1+e2+e3+e4-e5)=207*SQRT(3266594)/372775-23*SQRT(5037200385)/7012200

"The second term is not a basic integral, but a function that can be evaluate~
d."

int5 SUB 2

-50082214766438400*AA(e1+e2+e3+e4-e6)/5073574307278031

"The function is:"

AA(e1+e2+e3+e4-e6)=AF([1,1,1,1,0,-1],4,6,an,bt,1,0)

AA(e1+e2+e3+e4-e6)=28*SQRT(3266594)/54405-14*SQRT(5037200385)/4805775

"The third term is not a basic integral, but a function that can be evaluated~
."

int5 SUB 3

3084876849714687749721480000*AA(e2+e3+e4)/98708359449876165832998833

"The function is:"

AA(e2+e3+e4)=AF([0,1,1,1,0,0],4,6,an,bt,1,0)

AA(e2+e3+e4)=45*SQRT(3266594)/47554-SQRT(5037200385)/103974

"The fourth term is a basic integral."

int5 SUB 4

-33649773397162691621392678784*IA(-e1)/866165854172663355184564759575

"The basic integral is:"

IA(-e1)

IA(-e1)

"The fifth term is a basic integral."

int5 SUB 5

-14335213855710841650792*IA(-e5)/86521330480155184055

"The basic integral is:"

IA(-e5)

IA(-e5)

"The sixth term is a basic integral."

int5 SUB 6

23131892295883231600208*IA(-e6)/93379135125452160555

"The basic integral is:"

IA(-e6)

IA(-e6)

"I(-e6) and I(-e5) are the same basic integral with different parameters."

IA(-e6)=SUBST(IA(-e5),[a5],[a6])

IA(-e6)=IA(-e5)

"The siventh term is a basic integral."

int5 SUB 7

-12099952509405760834240952*IA(0)/320597342132401837783107

"The basic integral is:"

IA(0)

IA(0)

"Find the basic integrals in a table of integrals."

IS([IA(-e1),IA(-e5),IA(-e5),IA(0)],4)

[I428(a,b,x,y),I427(a,b,x,y),I427(a,b,x,y),I414(a,b,x,y)]

"Evaluate the basic integrals."

int6:=APPROX([I428(an,bt,1,0),I427(an,bt,1,0),I427([an SUB 1,an SUB 2,an SUB ~
3,an SUB 4,an SUB 6],bt,1,0),I414(an,bt,1,0)])

int6:=[882648970401909281320937426536171/770139438851967068305953537147127,75~
8959740946133335799407902352113/1558210957233881835349441486493482,5994369562~
10769012253723964997082/1363393030911334614424585861489257,103112358868101151~
4824204915500365/1062983952718432652285585411008329]

"Put the values of the A functions in an array."

int7:=[[e1+e2+e3+e4-e5,AF([1,1,1,1,-1,0],4,6,an,bt,1,0)],[e1+e2+e3+e4-e6,AF([~
1,1,1,1,0,-1],4,6,an,bt,1,0)],[e2+e3+e4,AF([0,1,1,1,0,0],4,6,an,bt,1,0)]]

int7:=[[e1+e2+e3+e4-e5,207*SQRT(3266594)/372775-23*SQRT(5037200385)/7012200],~
[e1+e2+e3+e4-e6,28*SQRT(3266594)/54405-14*SQRT(5037200385)/4805775],[e2+e3+e4~
,45*SQRT(3266594)/47554-SQRT(5037200385)/103974]]

"Put the values of the basic integrals in an array."

int8:=[[-e1,int6 SUB 1],[-e5,int6 SUB 2],[-e6,int6 SUB 3],[0,int6 SUB 4]]

int8:=[[-e1,882648970401909281320937426536171/7701394388519670683059535371471~
27],[-e5,758959740946133335799407902352113/1558210957233881835349441486493482~
],[-e6,599436956210769012253723964997082/1363393030911334614424585861489257],~
[0,1031123588681011514824204915500365/1062983952718432652285585411008329]]

"Substitute the values of the basic integrals and A functions values."

AA(v):=(SELECT(k SUB 1=v,k,int7)) SUB 1 SUB 2

IA(v):=(SELECT(k SUB 1=v,k,int8)) SUB 1 SUB 2

"Evaluate the integral."

int9:=APPROX(int4)

int9:=471935107228402336905440604717785/8664424068785748153380169455473276

"Compare results to numerical integration."

APPROX(int9-int2)

-3.056665831142513982007817197092306333057452041262191816753492316*10^(-15)

"The results match."

"Compare results to MAPLE."

"Generate the rational expression."

int10:=1/(SQRT((t+a1)*(t+a2)*(t+a3)*(t+a4))*(t+a5)^2*(t+a6)^2)

int10:=1/((a5+t)^2*(a6+t)^2*SQRT((a1+t)*(a2+t)*(a3+t)*(a4+t)))

"Factor the polynomial so MAPLE can integrate it easier."

int11:=FACTOR(SUBST(int10,[a1,a2,a3,a4,a5,a6],an))

int11:=1455723360*SQRT(465)/((28*t+53)^2*(69*t+116)^2*SQRT((65*t+29)*(75*t+43~
)*(117*t+73)*(124*t+119)))

"Digits:=64;"

"evalf(int(1455723360*sqrt(465)/((28*t+53)^2*(69*t+116)^2*sqrt((65*t+29)*(75*~
t+43)*(117*t+73)*(124*t+119))),t=0..1));"

".5446814508174694901018489489557609990731103509573270910308616e-1"

int12:=0.05446814508174694901018489489557609990731103509573270910308616

int12:=680851813521836862627311186194701248841387938696658863788577/125000000~
00000000000000000000000000000000000000000000000000000

APPROX(int9-int12)

8.850806991276468828638079999048891842430218924724353443574233445*10^(-62)

"The results match to 61 digits."