Riemannův integrál

Tomáš Kalvoda, KAM FIT ČVUT, 2013

V tomto notebooku si ukážeme na příkladu funkce f(x)=2+(1/5)(x^2)-x na intervalu ⟨0,6⟩. Zvídavý čtenář může snadno experimentovat s jinými volbami.

f[x_]=2+(1/5)(x^2)-x;
a=0;
b=6;

Graf této funkce je uveden na následujícím obrázku. Riemannův (určitý) integrál funkce f nad intervalem ⟨a,b⟩ pak představuje obsah plochy pod grafem funkce f.

Plot[f[x],{x,a,b},
PlotStyle->{Thick,Blue},
PlotRange->{-.5,3.5},
Filling->Axis,AxesLabel->{Style["x",14,Italic],Style["y",14,Italic]}
]

Pozor, plocha je chápána jako orientovaná. Tedy nad osou x s kladným znamínkem a pod osou x se záporným znaménkem. Například integrál z funkce sin na intervalu ⟨-π/2,π/2⟩ je nulový. Podívejte se na obrázek a výsledek integrace:

Plot[Sin[x],{x,-π/2,π/2},
PlotStyle->{Thick,Blue},
PlotRange->{-1.2,1.2},
Filling->Axis,AxesLabel->{Style["x",14,Italic],Style["y",14,Italic]}
]
Integrate[Sin[x],{x,-π/2,π/2}]

Dolní a horní součty

Uvažme rozdělení intervalu ⟨a,b⟩ na n∈ℕ ekvidistantních (stejně dlouhých) dílků. Každý dílek má tedy délku Δ=((b-a)/2). Rozdělení značínem řeckým písmenkem σ.

n=10;
delta=(b-a)/n;
sigma=Table[a+(i-1)delta,{i,1,n+1}];

Příslušným horním a dolním součtem jsou čísla:

minima=Table[
NMinValue[{f[x],sigma[[i]]<x<sigma[[i+1]]},x],
{i,1,Length[sigma]-1}
];
dolni=Sum[
m*delta,
{m,minima}
];
maxima=Table[
NMaxValue[{f[x],sigma[[i]]<x<sigma[[i+1]]},x],
{i,1,Length[sigma]-1}
];
horni=Sum[
m*delta,
{m,maxima}
];

A jejich grafické reprezentace.

Print["Horní součet: ",horni]
Plot[f[x],{x,a,b},
PlotStyle->{Thick,Blue},
PlotRange->{-.5,3.5},
Prolog->{Lighter[Blue],Opacity[.5],EdgeForm[Dashed]}~Join~Table[
Rectangle[{a+(i-1)*delta,0},{a+i*delta,maxima[[i]]}],{i,1,Length[maxima]}
],
AxesLabel->{Style["x",14,Italic],Style["y",14,Italic]}
]//Print
Print["Dolní součet: ",dolni]
Plot[f[x],{x,a,b},
PlotStyle->{Thick,Blue},
PlotRange->{-.5,3.5},
Prolog->{Lighter[Blue],Opacity[.5],EdgeForm[Dashed]}~Join~Table[
Rectangle[{a+(i-1)*delta,0},{a+i*delta,minima[[i]]}],{i,1,Length[minima]}
],
AxesLabel->{Style["x",14,Italic],Style["y",14,Italic]}
]//Print

Na závěr se zbavme nepotřebných proměnných (níže by mohlo dojít ke kolizi).

Clear[n,delta,sigma,minima,dolni,maxima,horni]

Na závěr ještě interaktivní manipulace (nezapomeňte pro jistotu na Clear).

Manipulate[
delta=(b-a)/n;
sigma=Table[a+(i-1)delta,{i,1,n+1}];
extremy=Table[
func[{f[x],sigma[[i]]<x<sigma[[i+1]]},x],
{i,1,Length[sigma]-1}
];
soucet=delta*Plus@@extremy;
Column[{
"Hodnota součtu: "~~ToString[soucet],
Plot[f[x],{x,a,b},
PlotStyle->{Thick,Blue},
PlotRange->{-.5,3.5},ImageSize->{500,400},
Prolog->{Lighter[Blue],Opacity[.5],EdgeForm[Dashed]}~Join~Table[
Rectangle[{a+(i-1)*delta,0},{a+i*delta,extremy[[i]]}],{i,1,Length[extremy]}
],
AxesLabel->{Style["x",14,Italic],Style["y",14,Italic]}
]},Alignment->Center
],
{{n,4,Style["n",Italic]},2,30,1,ContinuousAction->False,Appearance->"Labeled"},
{{func,NMinValue,""},{NMinValue->"Dolní součty",NMaxValue->"Horní součty"}}
]

Primitivní funkcí k funkci f je funkce F, kterou snadno zjitíme (i ručně).

F[x_]=Integrate[f[x],x]

Tudíž plocha vypočítávaná na obrázku výše je přesně:

N[F[b]-F[a],10]

Integrální součet

Jak víme z přednášek, přibližnou hodnotu určitého integrálu nemusíme počítat pomocí maxim a minim, která je složitá hledat, ale můžeme náhodně volit hodnotu funkce v bodě dělícího intervalu. Pro jednoduchost volme třeba jeho střed.

Manipulate[
delta=(b-a)/n;
sigma=Table[a+(i-1)delta,{i,1,n+1}];
hodnoty=Table[
If[vyber=="prumery",
f[(sigma[[i+1]]+sigma[[i]])/2],
f[RandomReal[{sigma[[i]],sigma[[i+1]]}]]
],
{i,1,Length[sigma]-1}
];
soucet=delta*Plus@@hodnoty;
Column[{
"Hodnota integrálního součtu: "~~ToString[N[soucet]],
Plot[f[x],{x,a,b},
PlotStyle->{Thick,Blue},
PlotRange->{-.5,3.5},ImageSize->{500,400},
Prolog->{Lighter[Blue],Opacity[.5],EdgeForm[Dashed]}~Join~Table[
Rectangle[{a+(i-1)*delta,0},{a+i*delta,hodnoty[[i]]}],{i,1,Length[hodnoty]}
],
AxesLabel->{Style["x",14,Italic],Style["y",14,Italic]}
]},Alignment->Center
],
{{n,4,Style["n",Italic]},2,30,1,ContinuousAction->False,Appearance->"Labeled"},
{{vyber,"prumery",""},{"nahodne"->"Náhodně","prumery"->"Průměry (středy)"}}
]

Lichoběžníkové pravidlo

V případě braní průměrné hodnoty ale nejde o nic jiného než o sčítání obsahů lichoběžníků. Pokud si představíte lichoběžník s podstavou Δ a svislými stranami dlouhými a a b, a<b, pak jeho obsah je

S=a·Δ+(1/2)(b-a)·Δ=Δ·(1/2)(a+b)

Graphics[{
Lighter[Blue],EdgeForm[Dashed],
Polygon[{{0,0},{1,0},{1,3},{0,2.5}}],
Black,
Text[Style["a",Italic,16],{-0.2,1.25}],
Text[Style["b",Italic,16],{1.2,1.5}],
Text[Style["Δ",Italic,16],{0.5,-.2}]
}
]

Manipulate[
delta=(b-a)/n;
sigma=Table[a+(i-1)delta,{i,1,n+1}];
hodnoty=Table[
f[(sigma[[i+1]]+sigma[[i]])/2],
{i,1,Length[sigma]-1}
];
soucet=delta*Plus@@hodnoty;
Column[{
"Hodnota integrálního součtu: "~~ToString[N[soucet]],
Plot[f[x],{x,a,b},
PlotStyle->{Thick,Blue},
PlotRange->{-.5,3.5},ImageSize->{500,400},
Prolog->{Lighter[Blue],Opacity[.5],EdgeForm[Dashed]}~Join~Table[
Polygon[{{a+(i-1)*delta,0},{a+i*delta,0},{a+i*delta,f[sigma[[i+1]]]},{a+(i-1)*delta,f[sigma[[i]]]}}],{i,1,Length[sigma]-1}
],
AxesLabel->{Style["x",14,Italic],Style["y",14,Italic]}
]},Alignment->Center
],
{{n,4,Style["n",Italic]},2,30,1,ContinuousAction->False,Appearance->"Labeled"}
]

Simpsonovo pravidlo

Na dělícím intervalu ⟨c,d⟩ nahraďme integrand kvadratickou funkcí (ne přímkou jako výše) procházející body (c,f(c)), (d,f(d)) a průměrem (m,f(m)), kde m=((c+d)/2). Následuje odvození vzorce.

inter[x_]=aa (x^2)+bb x+cc;
solution=Solve[
{inter[c]==ff[c],
inter[(c+d)/2]==ff[(c+d)/2],
inter[d]==ff[d]
},
{aa,bb,cc}
]//First

Integrujme tuto funkce na ⟨c,d⟩, dostáváme:

Integrate[inter[x]/.solution,{x,c,d}]

Dostáváme tak Simpsonovo pravidlo, které spočívá v nahrazení integrálu nad intervalem ⟨c,d⟩ integrálem z tohoto kvadratického interpolantu. Konkrétně

(∫_c^d)f(x)dx≈(1/6)(d-c)(f(c)+f(d)+4f(((c+d)/2))).

f[x_]=2+Sin[4x];
goal=Integrate[f[x],{x,a,b}];
Manipulate[
delta=(b-a)/n;
sigma=Table[a+(i-1)delta,{i,1,n+1}];
simpson=Table[
(1/6)delta(f[sigma[[i]]]+f[sigma[[i+1]]]+4f[(sigma[[i]]+sigma[[i+1]])/2]),
{i,1,Length[sigma]-1}
];
soucet=Plus@@simpson;
Column[{
"Hodnota součtu podle Simpsona: "~~ToString[N[soucet]],
"Skutečná hodnota: "~~ToString[N[goal]],
Plot[f[x],{x,a,b},
PlotStyle->{Thick,Blue},
PlotRange->{-.5,3.5},ImageSize->{500,400},
AxesLabel->{Style["x",14,Italic],Style["y",14,Italic]},
Prolog->{Gray,Dashed}~Join~Table[
Line[{{sigma[[i]],0},{sigma[[i]],f[sigma[[i]]]}}],
{i,1,Length[sigma]}]
]},Alignment->Center
],
{{n,4,Style["n",Italic]},2,30,1,ContinuousAction->False,Appearance->"Labeled"}
]

Není jednoduché napočítávat kvadratické interpolace, proto zobrazujeme jenom šířky sloupců pro ilustraci.