Numerikus Integrálás

Alapfeladat: (∫_a^b)f[x]dx Riemann integrál közelítése.

Alapötletek: 
1) (∫_a^b)f[x]dx≈(∫_a^b)L[x]dx, ahol L f interpolációs polinomja.

2) Ekvidisztáns alappontrendszer a≤(Subscript[x, 0])<…<(Subscript[x, n-1])≤b, (Subscript[x, k])-(Subscript[x, k-1])=h

Emlékeztető: (∫_a^b)f[x]dx≈(∫_a^b)L[x]dx=(∫_a^b)(∑_(i=0)^(n-1))f[(Subscript[x, i])](Subscript[l, i])[x]dx=(∑_(i=0)^(n-1))f[(Subscript[x, i])](∫_a^b)(Subscript[l, i])[x]dx=(∑_(i=0)^(n-1))f[(Subscript[x, i])](Subscript[w, i])

Kvadratúraformula:

(Subscript[Q, n])[f]=(∑_(i=0)^(n-1))f[(Subscript[x, i])](Subscript[w, i])

Interpolációs kvadratúraformula esetén a súlyok a Lagrange bázispolinomok integráljai.

Lagrange előállítása: bázispolinomok lineáris kombinációja

LagrBase[j_,x_List,var_]:=((∏_(k=1)^j)((var-x[[k]])/(x[[j+1]]-x[[k]])))((∏_(k=j+2)^(Length[x]))((var-x[[k]])/(x[[j+1]]-x[[k]])))

LagrInterp[x_List,y_List,var_]:=(∑_(j=0)^(Length[x]-1))y[[j+1]]LagrBase[j,x,var]

Trapéz (formula)

A formula levezetése, pontossági rendje

LagrInterp[{a,b},{1,0},x]

LagrBase[0,{a,b},x]

InterpolatingPolynomial[{{a,1},{b,0}},x]

Together[%]

Integrate[LagrBase[0,{a,b},x],x]

(%/.x->b)-(%/.x->a)//Simplify

Integrate[LagrBase[0,{a,b},x],{x,a,b}]//Simplify

Integrate[LagrBase[1,{a,b},x],{x,a,b}]//Simplify

Formula:

Q2=(b-a)/2(fa+fb)

Rend:

Integrate[1,{x,a,b}]

Q2/.{fa->1,fb->1}

Integrate[x,{x,a,b}]//Factor

Q2/.{fa->a,fb->b}

Integrate[(x^2),{x,a,b}]//Factor

Q2/.{fa->(a^2),fb->(b^2)}

Integrate[Binomial[t,2],{t,0,1}]

Integrate[x^2,{x,1,2}]

g1=Plot[x^2,{x ,1,2},FillingStyle->Red,Filling->Axis]

N[{7/3,5/2}]

g2=Plot[3x-2,{x ,1,2},FillingStyle->Yellow,Filling->Axis]

Show[g2,g1]

(∫_a^b)LagrInterp[{a,b},{fa,fb},x]\[DifferentialD]x

Integrate[x^2,{x,1,2}]

LagrInterp[{1,2},{1,4},x]//Expand

Integrate[%,{x,1,2}]//N

Plot[{x^2,3x-2},{x,1,2},AxesOrigin->{1,0},PlotRange->{0,9/2},FillingStyle->Automatic,Filling->Axis,AspectRatio->1.5,ImageSize->{300,300}]

Simpson formula

LagrInterp[{a,(a+b)/2,b},{fa,fm,fb},x]

(∫_a^b)LagrInterp[{a,(a+b)/2,b},{fa,fm,fb},x]\[DifferentialD]x

Q3=%;

LagrInterp[{-1,0,1},{1,0,1},x]//Expand

Plot[(x^2),{x,-2,2},Filling->Automatic]

Végezzük el az alábbi lépésket, ha az(∫_1^2) (1/x)dx határozott integrált közelítjük!

f[x_]:=1/x;
{a1,b1}={1,2};

Q2

Q2/.{a->a1,b->b1,fa->f[a1],fb->f[b1]}

Q3

Q3/.{a->a1,b->b1,fa->f[a1],fb->f[b1],fm->f[(a1+b1)/2]}

Végezzük el az alábbi lépésket, ha  a (∫_0^2) (e^(-(x^2)))dx határozott integrált közelítjük!

g[x_]:=E^(-x^2);
{a2,b2}={0,2};

(∫_0^2) (e^(-(x^2)))dx  határozott integrált közelítése:

g1a=Plot[g[x],{x,a2,b2},Filling->Axis]

g1b=Plot[LagrInterp[{a2,(a2+b2)/2,b2},{g[a2],g[((a2+b2)/2)],g[b2]},x],{x,a2,b2},Filling->Axis,FillingStyle->Yellow]

Show[g1a,g1b]

g1b=Plot[{g[x],LagrInterp[{a2,(a2+b2)/2,b2},{g[a2],g[((a2+b2)/2)],g[b2]},x]},{x,a2,b2},Filling->Axis,FillingStyle->Automatic]

A pontos érték (Newton-Leibniz szabály)

(∫_a2^b2)g[x]\[DifferentialD]x

N[%]

Beépített numerikus integrálás:

NIntegrate[g[x],{x,a2,b2}]

Az interpolációs polinom integrálása

(∫_a2^b2)LagrInterp[{a2,(a2+b2)/2,b2},{g[a2],g[((a2+b2)/2)],g[b2]},x]\[DifferentialD]x

N[%]

A Simpson formula direkt alkalmazása

((b2-a2)/6)(g[a2]+4 g[((a2+b2)/2)]+g[b1])

N[%]

Simpson szabály

Alapötlet: Bontsuk a “hosszú” [a,b] intervallumot részintervallumokra és minden részintervallumra alkallmazzuk a Simpson formulát!
Ekkor I~(∑_j)(Subscript[Q, s])[f,(Subscript[I, j])].
Simpson súlyfüggvény:

Clear[WS];
WS[0,_]=1;
WS[m_,m_]=1;
WS[n_?EvenQ,_]:=2;
WS[n_?OddQ,_]:=4;

Sum[WS[i,4] F[i],{i,0,4}]

Sum[WS[i,10] F[i],{i,0,10}]

(∫_1^2) (1/x)dx határozott integrált közelítése (f primitiválható)

f[x]

m=10, h=1/10 alappontok:  (Subscript[x, k])=1+k/10; (k=0,1,2,…,10)

m=100;h=(1/100);

(∑_(k=0)^m)F[1+k/100]WS[k,m]

(h/3)(∑_(k=0)^m)f[1+k/100]WS[k,m]

N[%]

N[Log[2]]

(∫_1^2)f[x]\[DifferentialD]x

Integrate[x^2,{x,1,2}]

g1=Plot[x^2,{x ,1,2},FillingStyle->Red,Filling->Axis]

N[{7/3,5/2}]

g2=Plot[3x-2,{x ,1,2},FillingStyle->Yellow,Filling->Axis]

Show[g2,g1]

Az m növelésével a hiba csökkenthető! 1/180((((b-a))^5)/(m^4)) (f^((IV)))[ξ]

Zárt Newton-Cotes formulák

Clear[h]

(Subscript[x, 0])=a,(Subscript[x, n-1])=b, h=((b-a)/n)

EDPoints[sp_,n_,h_]:=Table[sp+(i)h,{i,0,n-1}]

FPoints[f_,n_]:=Table[ToExpression[StringJoin[ToString[f],ToString[i]]],{i,0,n-1}]

EDPoints[a,5,h]

FPoints[f,5]

LagrInterp[EDPoints[a,2,h],FPoints[f,2],x]

Factor[Integrate[LagrInterp[EDPoints[a,2,h],FPoints[f,2],x],{x,a,a+h}]]

Figyelem: h minden sorban más!

TableForm[Table[Factor[Integrate[LagrInterp[EDPoints[a,n,h],FPoints[f,n],x],{x,a,a+(n -1)h}]],{n,2,10}]]

Itt a Simpson formula (második a sorban):

SF=Factor[Integrate[LagrInterp[EDPoints[a,3,h],FPoints[f,3],x],{x,a,a+h,a+2h}]]/.{h->(b-a)/2}

Súlyok (weights):

Map[Coefficient[SF,#]&,{f0,f1,f2}]

Nyílt  Newton-Cotes formulák

Itt a végpontok nem tartoznak az alappontokhoz.

TableForm[Table[Factor[Integrate[LagrInterp[Drop[EDPoints[a,n-1,h],1],Drop[FPoints[f,n-1],1],x],{x,a,a+(n -1)h}]],{n,3,10}]]

Hibabecslés és pontossági rend

Pontossági rend: Ellenőrizzük, hogy a Simpson formula esetén r=3!

polylist={1,x,(x^2),(x^3),(x^4)};

(∫_a^b)polylist\[DifferentialD]x

Table[((b-a)/6) ((a^n)+4(((((a+b)/2)))^n)+(b^n)),{n,0,4}]//Expand

A Simpson formula alkalmazásával közelítjük a (∫_0^1)(e^(-(x^2)))dx integrált. Becsüljük meg a hibát!

(((b-a))^5)(f^( (4)))[ξ]/2880

32 90

t=Table[Factor[D[(E^(-(x^2))),{x,n}]],{n,0,4}]

Nagyon durva becslése a negyedik deriváltnak (76), intervallum hossza 1.

76/(90 32)//N

Plot[Abs[Last[t]],{x,0,1}]

NIntegrate[E^(-x^2),{x,0,1}]

1/6(1+4 E^(-1/4)+1/E)//N

Általánosított kvadratúra probléma: (∫_a^b)f[x]ω [x]dx
Most csak azzal foglakozunk, amikor ω=1, [a,b]=[-1,1].
Nem ekvidisztáns alappontrendszer, n pont esetén [-1,1]-en minden(!) (2n-1)-ed fokú polinomra pontos;

Speciálisan:
	n=2, 3-adfokú polinomig pontos
	n=3, 5-ödfokú polinomig pontos

Hogyan kapjuk a pontokat és a súlyokat? (∫_a^b)f[x]dx=(∑_(i=0)^(n-1))f[(Subscript[x, i])](Subscript[w, i])

Alappontok: Ortogonális polinomrendszer elemeinek gyökei
Súlyok: n=2 esetén {1,1} (nem kell tudni általánosan)
Ha <f|g>=(∫_(-1)^1)f[x]g[x]dx, akkor az  OPR elemei az ún. Legendre polinomok.

Speciálisan,  n=2 esetén a Legendre polinom (1/2) (3 (x^2)-1).

Kísérlet a súlyok és az osztópontok meghatározására a pontossági rendre vonatkozó feltételből.
n=2: osztópontok x[0], x[1]
n=3: osztópontok x[0], x[1], x[2]

Clear[W,X,w,x];

W=Table[w[j],{j,0,5}];
X=Table[x[j],{j,0,5}];

sop[expr_,var_,n_]:=Sum[W[[j]] expr/.x->X[[j]],{j,n}]

t1=Table[sop[x^n,x,2],{n,0,3}]

t2=Table[Integrate[x^n,{x,-1,1}],{n,0,3}]

Thread[Equal[t1,t2]]

Solve[Thread[Equal[t1,t2]]]

t1=Table[sop[x^n,x,3],{n,0,5}]

t2=Table[Integrate[x^n,{x,-1,1}],{n,0,5}]

Solve[Thread[Equal[t1,t2]]][[4]]

Megjegyzés: Az osztópontok permutációtól eltekintve egyértelmû.

Érdekesség n=4 esetén:

t1=Table[sop[x^n,x,4],{n,0,7}];

t2=Table[Integrate[x^n,{x,-1,1}],{n,0,7}];

GroebnerBasis[Thread[Plus[t1,-t2]],Join[Take[X,4],Take[W,4]]]

Solve[%[[1]]==0]

N[%]

Speciális eset: a=-1,b=1, ρ=1, n=3

q[x_]:=c0+c1 x+c2 (x^2)+c3 (x^3);

q ortogonális  az 1,x,(x^2) polinomokra:

{(∫_(-1)^1) q[x] 1\[DifferentialD]x==0,(∫_(-1)^1) q[x] x\[DifferentialD]x==0,(∫_(-1)^1) q[x] (x^2)\[DifferentialD]x==0}

Solve[%]

q[x]/.{{c0->0,c2->0,c1->-((3 c3)/5)}}/.{c3->5/2}

Ez a harmadrendű Legendre polinom:

LegendreP[3,x]

?LegendreP

LegendreP[2,x]

A Legendre polinomok ortogonális rendszert alkotnak.
Az ortogonalitás ellenőrzése: <(Subscript[L, 1])|(Subscript[L, 2])> =0

Integrate[LegendreP[1,x] LegendreP[2,P],{x,-1,1}]

TableForm[Table[{i,LegendreP[i,x]},{i,0,10}]]

Ortogonalitás:

TableForm[Table[Integrate[LegendreP[i,x] LegendreP[j,x],{x,-1,1}],{i,0,5},{j,0,5}]]

Mintapontok: a Legendre polinomok gyökei

Solve[LegendreP[3,x]==0]

r3=Sort[x/.Solve[LegendreP[3,x]==0],Less]

A súlyok meghatározása a pontossági rend alapján.

plist={1,x,(x^2),(x^3),(x^4),(x^5)};

Table[(∑_i^3)W[[i]] (plist[[j]]/. x->r3[[i]])==(∫_(-1)^1)plist[[j]]\[DifferentialD]x,{j,3}]

Solve[%]

Itt vannak a súlyok általában  (2/(1-((Subscript[x, i])^2)((P ' [(Subscript[x, i])])^2))):

Weight[i_,n_,var_]:=Module[{xi=Sort[var/. Solve[LegendreP[n,var]==0],Less][[i+1]]},(2/((1-(xi^2)) ((((∂_var)LegendreP[n,var]/. var->xi))^2)))]

Simplify[Weight[2,3,x]]

Simplify[Table[Weight[i,n,x],{n,5},{i,0,n-1}]]

Speciális eset: a=-1,b=1, ρ=1, n=3

Map[sop[#,x,3]&,Table[LegendreP[j,x],{j,0,2}]]/.{x[0]->-Sqrt[3/5],x[1]->0,x[2]->Sqrt[3/5]}

Thread[Equal[%,{2,0,0}]]

Solve[%]

Hasonlítsuk össze a trapézszabályt, az n=2-höz tartozó Gauss kvadratúrát és az integrál pontos értékét, ha , f(x)=(x^3)+(x^2)-1

f1[x_]:=(x^3)+(x^2)-1;

(∫_(-1)^1)f1[x]\[DifferentialD]x

p=Integrate[f1[x],x]

(p/.x->1)-(p/.x->-1)

(1/2)(f1[-1]+f1[1])

Nem pontos!

Solve[LegendreP[2,x]==0]

1 f1[-1/Sqrt[3]]+1 f1[1/Sqrt[3]]

A Gauss kvadratúra pontos értéket ad, a formula pontossági rendje 3, vagyis harmadfokú polinomokra a hiba 0!

f2[x_]:=(x^4);

f2[2]

1 f2[-1/Sqrt[3]]+1 f2[1/Sqrt[3]]

(∫_(-1)^1)f2[x]\[DifferentialD]x