Izbira koraka

<<Numerična matematika

Trapezno pravilo s kontrolo napake

Približek za napako T(h) − I

Približek za napako T(h / 2) − I

Če želimo izračunati približek za napako, najprej izračunamo:

in nato še

Približek za napako:

Če je približek za napako po absolutni vrednosti manjši od predpisane natančnosti , je T((b − a) / 2) dovolj dober približek za I, sicer korak h razpolovimo in izračun ponovimo.

Algoritem: Trapezno pravilo s kontrolo napake

Naj bo f zvezna funkcija na [a,b], N neko naravno število in . Naslednji algoritem izračuna s pomočjo trapezne formule približek T, ki se od vrednosti določenega integrala razlikuje za manj kot ali pa se po N razpolovitvah dolžine koraka konča brez rezultata (T=NaN).

e=2*eps
m=0
h=b-a
T=h*(f(a)+f(b))/2

while (m<N) & < (abs(e) > eps)
     m=m+1
     h=h/2
     k=2^(m-1)
     s=0
     for i=1:k
          s=s+f(a+(2*i-1)*h)
     end
     e=s*h-T/2
     T=T+e
end

if abs(e) > eps
     T=NaN
end

Trapezno pravilo s adaptivno izbiro koraka

Vse dosedanje metode so uporabljale delitev na enake podintervale. Včasih je bolje celoten interval razdeliti na različno dolge intervale, katerih dolžino izbiramo glede na lokalno obnašanje funkcije (tisti deli intervala, kjer je f(x) manjši, vzamemo daljše podintervale, kjer je pa večji pa manjše).

Če naj bo napaka na celotnem intervalu [a,b] manjša od , je simselna zahteva, da naj bo napaka na podintervalu dolžine h_i manjša od . Delni rezultat T(h_i / 2) torej sprejmemo, če je

sicer pa moramo vzeti manjši podinterval. Na osnovi ocene napake tudi lahko določimo optimalno dolžino naslednjega podintervala:

kjer smo s σ označili "varnostni koeficient", ki ga navadno izberemo malo manj kot 1 (npr. 0,9).

Algoritem: Adaptivno trapezno pravilo

Naj bo f zvezna funkcija na intervalu [a,b] in . Naslednji algoritem izračuna s pomočjo trapezne formule približek I, ki se od vrednosti določenega integrala razlikuje manj kot .

sigma=0.9
x=a
I=0
h=b-a
f_x=f(a)
f_xh=f(b)

while x<b
     f_xh2=f(x+h/2)
     T_1=h*(f_x+f_xh)/2          % Najprej z osnovnim korakom
     T_2=T_1/2 + h*f_xh/2/2      % Ponovno, s polovičnim korakom
     if (abs((T_1-T_2)/3) < h*eps/(b-a)
          x=x+h                  % rezultat sprejet
          I=I+T_2                % prištejemo delni rezultat
          f_x=f_xh
          h=sigma*sqrt(3*h^3*eps/abs((T_1-T_2)/(b-a))          % novi korak
          
          if x+h > b             % ali smo že blizu konca
               h=b-x
          end

          f_xh = f(x+h)
     else                        % rezultat zavrnjen
          h=h/2                  % razpolovimo korak
          f_xh=f_xh2
     end
end