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

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)
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

1 Like

DNF
January 16, 2023, 8:49am
5
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; ...)`

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.

DNF
January 16, 2023, 9:07am
7
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

DNF
January 16, 2023, 9:09am
8
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.