Make the substitution #u=x^2#. Notice #"d"x = 1/(2x) "d"u#.
Then,
#int x sqrt(x^4+a^4) "d"x = int x sqrt(x^4+a^4) 1/(2x) "d"u#
#int x sqrt(x^4+a^4) "d"x = 1/2 int sqrt(u^2+a^4) "d"u#
This leaves us with an integral which is really suited to be solved with a hyperbolic substitution. Substitute #u=a^2 sinh(t)#. Then,
#("d"u)/("d"t) a^2cosh(t)#,
#"d"u = a^2cosh(t) "d"t#.
Then, the integral is transformed to,
#1/2 int sqrt(a^4sinh^2(t) + a^4) a^2cosh(t) "d"t#
#a^2/2 int sqrt(a^4) sqrt(sinh^2(t) + 1) cosh(t) "d"t#.
Using the identity #cosh^2(x)-sinh^2(x)=1# gives #cosh(x)=sqrt(1+sinh^2(x))#.
Then,
#a^4/2 int cosh^2(t) "d"t#.
The identity #cosh(2x)=2cosh^2(x)-1# gives #cosh^2(x) = 1/2(cosh(2x)+1)#.
Then,
#a^4/2 int cosh^2(t) "d"t = a^4/4 int cosh(2t) + 1 "d"t#,
#a^4/4 int cosh(2t) +1 "d"t = a^4/8 sinh(2t) + a^4/4t + C#,
#a^4/4 int cosh(2t) +1 "d"t = a^4/4 sinh(t)cosh(t) + a^4/4t + C#.
As #u=a^2sinh(t)#, #t="arsinh"(u/a^2)#.
Then,
#a^4/4sinh(t)cosh(t)+a^4/4 t = a^4/4sinh(t)sqrt(sinh^2(t)+1) + a^4/4t#,
#a^4/4sinh(t)cosh(t) +a^4/4 t = a^4/4 sinh("arsinh"(u/a^2))sqrt(sinh("arsinh"(u/a^2))+1) + a^4/4"arsinh"(u/a^2)#,
#a^4/4sinh(t)cosh(t) +a^4/4 t = a^2/4 u sqrt(u^2/a^4+1)+a^4/4"arsinh"(u/a^2)#,
#a^4/4sinh(t)cosh(t) +a^4/4 t = u/(4) sqrt(u^2+a^4)+a^4/4"arsinh"(u/a^2)#.
The logarithmic definition of the inverse #sinh# function is #"arsinh"(x) = ln(x+sqrt(x^2+1))#.
#"arsinh"(u/a^2) = ln(u/a^2 + sqrt(u^2/a^4+1))#,
#"arsinh"(u/a^2) = ln((u+sqrt(u^2+a^4))/a^2)#,
#"arsinh"(u/a^2) = ln(u+sqrt(u^2+a^4))-ln(a^2)#.
We conclude that,
#a^4/4sinh(t)cosh(t) +a^4/4 t = u/(4) sqrt(u^2+a^4)+a^4/4ln(u+sqrt(u^2+a^4))-a^4/4ln(a^2)#.
Then, we need to reverse the final substitution, #u=x^2#.
#u/(4) sqrt(u^2+a^4)+a^4/4ln(u+sqrt(u^2+a^4))+a^4/4ln(a^2) = x^2/(4) sqrt(x^4+a^4)+a^4/4ln(x^2+sqrt(x^4+a^4))-a^4/4ln(a^2)#.
We conclude that,
#int x sqrt(x^4+a^4) "d"x = x^2/(4) sqrt(x^4+a^4)+a^4/4ln(x^2+sqrt(x^4+a^4))-a^4/4ln(a^2) + C#.
The #-a^3/4ln(a^2)# term can be absorbed by the constant #C# as it is just a constant.
We conclude,
#int x sqrt(x^4+a^4) "d"x = x^2/(4) sqrt(x^4+a^4)+a^4/4ln(x^2+sqrt(x^4+a^4))+C#.
