Hi all,
I have tried QuadGK, but the result is not the same as the solution manual.
Is there a problem in between?
this is my code:
using SymPy, QuadGK
x = symbols("x")
f(x) = (x^6 + 2)/(8x^2)
# simplify(diff(f(x),x)
g(x) = sqrt(1 + (x^6 - 2)/(2x^3))
h(x) = f(x)*g(x)
Area(x) = 2pi*h(x)
d = quadgk(x -> Area(x), 1, 3, rtol=1e-10)
println("c = ", d)
No its okay,
I forgot to power the derivative. Now the problem is solved.
But the result still differ a lot!
325.4437377507612 from QuadGK
If you use
g(x) = sqrt(1 + ((x^6 - 1)/(2x^3))^2)
the result is
c = (326.9195614457822, 1.2235672031124523e-8)
which compares very favorably to the exact answer
julia> 8429π/81
326.9195614457823
This is lovely:
using SymPy, QuadGK
x = symbols("x")
f(x) = (x^6 + 2)/(8x^2)
# simplify(diff(f(x),x))
# we cannot use diff inside the fd(x) since QuadGK can't comprehend diff
g(x) = sqrt(1 + ((x^6 - 1)/(2x^3))^2)
#g(x) = f(x)*√(1+(fd(x))^2)
h(x) = f(x)*g(x)
Area(x) = 2pi*h(x)
d = quadgk(x -> Area(x), 1, 3, rtol=1e-10)
dq = 2π*quadgk(h, 1, 3)[1]
println("(Area, error) = ", d)
println("Area with QuadGK = ", dq)
(Area, error) = (326.9195614457822, 1.2235672031124523e-8)
Area with QuadGK = 326.91956144578495
I’m not sure if you are aware, but the x = symbols("x") doesn’t have any effect here. The x inside the definitions of f(x), g(x), etc. are completely unrelated to the symbolic x that you create first. Its presence is a bit confusing and it doesn’t seem to belong in a Minimal Working Example.
One thing you could do, is to use automatic differentiation:
using ForwardDiff: derivative
f(x) = (x^6 + 2)/(8x^2)
A(x) = 2π * f(x) * sqrt(1 + derivative(f, x)^2)
julia> quadgk(A, 1, 3; rtol=1e-10)
(326.9195614457822, 1.2235672031124523e-8)
The autodiff version is only slightly slower:
julia> @btime quadgk(Area, 1, 3, rtol=1e-10)
897.619 ns (2 allocations: 368 bytes)
(326.9195614457822, 1.2235672031124523e-8)
julia> @btime quadgk(A, 1, 3, rtol=1e-10)
1.150 μs (2 allocations: 368 bytes)
(326.9195614457822, 1.2235672031124523e-8)
(It’s quite incredible how little performance overhead there is, actually.)
BTW, there’s no need to write x -> Area(x), this is just the same thing as Area. So in your code you would write: quadgk(Area, 1, 3; ...)
Well, first I am not aware, I thought using
x = symbols("x")
is needed to calculate the integral.
I never use Automatic Differentiation or try ForwardDiff, I will test and try it now.
In fact, if you re-write your f(x) to the equivalent form (taken from your posted equations)
f(x) = x^4/8 + 1/4x^2 # this is the same
A(x) = 2π * f(x) * sqrt(1 + derivative(f, x)^2)
then it’s even faster:
julia> @btime quadgk(A, 1, 3; rtol=1e-10)
696.528 ns (2 allocations: 368 bytes)
(326.91956144578234, 1.2235728874543383e-8)
The answer is slightly different from before, but still as accurate as the previous version.
No, it does nothing here, and is completely unrelated to quadgk, which does purely numerical integration.
It takes time for me to read and comprehend all your responses, thanks anyway.