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!

1 Like

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
``````
2 Likes

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)

println("(Area, error) = ", d)

println("Area with QuadGK = ", dq)
``````

(Area, error) = (326.9195614457822, 1.2235672031124523e-8)

1 Like

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)

(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; ...)`

3 Likes

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.

2 Likes

No, it does nothing here, and is completely unrelated to quadgk, which does purely numerical integration.

1 Like

It takes time for me to read and comprehend all your responses, thanks anyway.