Večrazsežni integrali

V poglavju o enojnih integralis smo spoznali, da je večina kvadraturnih formul preprosta utežena vsota

\[\int_a^b f(x)dx \approx \sum w_k f(x_k), \]

kjer so uteži $w_k$ in vozlišča $x_k$ izbrana tako, da je formula točna za polinome čim višjega reda.

Pri večkratnih integralih se stvari nekoliko zakomplicirajo, a v bistvu ostanejo enake. Kvadrature so tudi za večkratne integrale večinoma navadne utežene vsote vrednosti v izbranih točkah na območju.

Dvojni integral in integral integrala

Oglejmo si najbolj enostaven primer, ko integriramo funkcijo na kocki $[a,b]^2$. Dvojni integral lahko zapišemo kot dva gnezdena enojna integrala ^[1]

\[\int\int_{[a,b]^2} f(x,y)dxdy = \int_a^b\left(\int_a^b f(x,y)dy\right)dx = \int_a^b\left(\int_a^b f(x,y)dx\right)dy \]

Najbolj enostavno je izpeljati kvadrature za večkratni integral, če za vsak od gnezdenih enojnih integralov uporabimo isto kvadraturno formulo

\[\begin{equation}\label{eq:quad} \int_a^b f(x) dx \approx \sum_{k=1}^n w_k f(x_k) \end{equation}\]

z danimi utežmi $w_1, w_2, \ldots w_n$ in vozlišči $x_1, x_2, \ldots x_n$. Če za zunanji integral uporabimo kvadrature \eqref{eq:quad}, dobimo

\[\int\int_{[a,b]^2} f(x,y)dxdy = \int_a^b\left(\int_a^b f(x,y)dy\right)dx = \sum w_i Fy(x_i),\]

kjer je funkcija $Fy(x)$ enaka integralu po $y$, za katero lahko zopet uporabimo kvadrature \eqref{eq:quad}

\[Fy(x) = \int_a^b f(x,y) dy \approx \sum w_j f(x, y_j)\]

Dvojni integral lahko tako približno izračunamo kot dvojno vsoto

\[\begin{equation}\label{eq:quad2d} \int\int_{[a,b]^2} f(x,y)dxdy \approx \sum_{i,j} w_i w_j f(x_i, y_j). \end{equation}\]

Kvadraturni formuli, ki jo dobimo na ta način, pravimo produktna formula.

Dodatna naloga

Izpelji formulo za ostanek pri izračunu dvojnega integrala na kvadratu $[a,b]^2$, če za kvadraturo \eqref{eq:quad} uporabiš sestavljeno Simpsonovo pravilo

\[\int_a^b f(x) dx \approx \frac{h}{3} \sum_{i=0}^{n-1} \left(f(x_{2i}) + 4f(x_{2i+1}) + f(x_{2i+2})\right) + R_f,\]

kjer je $h=\frac{b-a}{2n}$, vozlišča $x_i = a +ih$ in ostanek je enak

\[R_f = -\frac{h^4}{180}(b-a)f^{(4)}(\xi),\]

za nek $\xi\in (a,b)$.

Povprečna razdalja med točkama na kvadratu $[0,1]^2$

Naloga

Izračunaj povprečno razdaljo med dvema točkama na kvadratu $[0,1]\times[0,1]$. Primerjaj različne metode.

Rešitev

Povprečna razdaljo lahko izračunamo s štirikratnem integralom

\[\bar{d} = \int_{[0,1]^4} \sqrt{(x_1-x_2)^2 + (y_1-y_2)^2} dx_1 dx_2 dy_1 dy_2.\]

Za izračun bom uporabili produktno kvadraturo s sestavljeno Simpsonovo formulo in metodo Monte Carlo.

using NumMat, LinearAlgebra
## povprečna razdalja med točkama v kocki [0,1]^2
f(x) = norm(x[1:2]-x[3:4]); # razdalja
x0 = LinRange(0, 1, 7); w = (x0[2]-x0[1])/3*[1 4 2 4 2 4 1];
I = ndquad(f, x0, w, 4)

0.5196600342442085

Poskusimo še z metodo Monte Carlo, kjer vozlišča izberemo slučajno enakomerno na izbranem območju.

function mcquad(fun, n, dim)
    Ef = 0
    for i=1:n
      Ef += fun(rand(dim))
    end
    return Ef/n
end
mcquad(f, 100000, 6)

0.5206771230765583

Primerjava različnih metod

Produktna kvadratura s sestavljeno Simpsonovo formulo

Poglejmo si, kako je s hitrostjo konvergence pri produktnih kvadraturah.

using Plots, Random
Random.seed!(1234)
ndquad_simpson(f, n, dim) = ndquad(f, LinRange(0, 1, 2n+1),
                            1/(6n)*vcat([1], repeat([4,2], n-1), [4, 1]), dim)
dim = 4
I = ndquad_simpson(f, 20, dim)
n = 3:15
napake_s = []
napake_mc = []
for nk in n
    push!(napake_s, ndquad_simpson(f, nk, dim) - I)
    push!(napake_mc, mcquad(f, (2nk+1)^4, dim) - I)
end
scatter((2n.+1).^4, abs.(napake_s), scale=:log10, label="simpson")
scatter!((2n.+1).^4, abs.(napake_mc), scale=:log10, label="Monte Carlo", title="Napake v odvisnosti od števila izračunov")

Z zbranimi podatki lahko določimo približni red simpsonove produktne metode za 4 kratne integrale

konst, red = hcat(ones(size(n)), log.((2n.+1).^4))\log.(abs.(napake_s))
println("Napaka produktne simpsonove formule pada z n^(", red, "), kjer je n število izračunov funkcijske vrednosti.")

Napaka produktne simpsonove formule pada z n^(-0.8871618142136483), kjer je n število izračunov funkcijske vrednosti.

Podobno lahko vsaj približno ocenimo hitrost konvergence za metodo Monte Carlo. Pri čemer se moramo zavedati, da je vrednost in tudi napaka odvisna od zaporedja slučajnih vozlišč, zato je ocena zgolj okvirna:

konst, red = hcat(ones(size(n)), log.((2n.+1).^4))\log.(abs.(napake_mc))
println("Napaka pri Monte Carlo pada približno z n^(", red, ") za izbrane vzorce.")

Napaka pri Monte Carlo pada približno z n^(-0.6398521320453916) za izbrane vzorce.

Koda

• NumMat.ndquad

I = ndquad(f, x0, utezi, d)

julia> f(x) = x[1] + x[2]; #f(x,y)=x+1;
julia> utezi = [1,1]; x0 = [0.5, 1.5]; #sestavljeno sredinsko pravilo
julia> ndquad(f, x0, utezi, 2)
8.0