Area of Surface of Revolution Integral too Hard to be Computed By JULIA' SymPy and Python' SymPy

General Usage
5

Let’s try numerical integration.

It takes three lines in Julia to make a trapezoidal integration function, out of which one line is end:

integrate(f, range) = let r=range, a=first(r), b=last(r), dx=step(r)
    (0.5f(a) + sum(f, r[2:end-1]) + 0.5f(b))dx
end

to check that it works:

julia> integrate(abs2, range(0,1,10^8))
0.3333333333333333

For the calculation of f^\prime(x) I’m going to use ForwardDiff, which uses dual numbers to automatically differentiate the function.

julia> using ForwardDiff
       Base.adjoint(f::Function) = Base.Fix1(ForwardDiff.derivative, f)

julia> f(x) = (x^6+2)/8x^2
f (generic function with 1 method)

julia> g(x) = f(x)*√(1+f'(x))
g (generic function with 1 method)

Ok let’s try integrating your function.

julia> 2π*integrate(g, range(1, 3, 10^8))
116.28129729349027

Let’s confirm against QuadGK, which uses magic to run much faster for the same precision (and even offers information about its numerical precision!)

julia> using QuadGK

julia> 2π*quadgk(g, 1, 3)[1]
116.2812972934902

Julia saves the day yet again (when does it not?).

3 Likes