Num. Math.—12

Gauss kvadratúra

Á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.

n=2, 3-adfokú
n=3, 5-ödfokú stb. 

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

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

2. Legendre polinom, pl. n=2:   1/2 (3x^2-1)

Határozatlan együtthatók módszere a formulák meghatározozására

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

2n-edfokú nemlineáris egyenletrendszer

TableForm[%]

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

Levezetett kvadratúraformula:

Q=f(-(1/(Sqrt[3])))+f(+(1/(Sqrt[3])))    ([a,b]=[-1,1])

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. Gröbner bázisok elmélete

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

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

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

Solve[%[[2,1]]==0]

N[w[3]/.%]

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

Solve[%[[2,1]]==0]

N[x[3]/.%]

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

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

AbsoluteTiming[First[GroebnerBasis[Thread[Plus[t1,-t2]],Join[Take[X,5],Take[W,5]]]]]

Solve[%[[2]]==0,w[4]]

N[%]

AbsoluteTiming[First[GroebnerBasis[Thread[Plus[t1,-t2]],Join[Take[W,5],Take[X,5]]]]]

Solve[%[[2]]==0,x[4]]

N[%]

Kísérlet (határozatlan együtthatók módszere ortogonális polinomok meghatározására): Speciális eset. a=-1,b=1, ρ=1, n=3

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

q ortogonális 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 köbös 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)}; w={w0,w1,w2};

Table[Sum[w[[i]](plist[[j]]/.x->r3[[i]]),{i,3}]==Integrate[plist[[j]],{x,-1,1}],{j,3}]

Solve[%]

Itt vannak a súlyok általában (2/(1-(Subscript[x, i])) (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)(D[LegendreP[n,var],var]/.var->xi)^2)]

Simplify[Weight[2,3,x]]

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

Kísérlet (osztópontok ismeretében súlyok) Speciális eset. a=-1,b=1, ρ=1, n=3

Clear[w]

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

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[%]

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

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

Manipulate[Plot[{0,LegendreP[n,x]},{x,-1,1},PlotRange->{-1,1},Filling->{1->{{2},{Red,Blue}}}],{n,1,7,1}]

Manipulate[Plot[{0,LegendreP[n,x]},{x,-1,1},PlotRange->{-1,1},Epilog->{Red,PointSize[.02],Point[{x,0}/.NSolve[LegendreP[n,x]==0,x]]}],{n,2,7,1}]

TT=Plot[Table[LegendreP[n,x],{n,1,11}],{x,-1,1},Epilog->{Red,PointSize[.02],Point[Flatten[Table[{x,0}/.NSolve[LegendreP[n,x]==0,x],{n,11}],1]]},PlotRange->{-1,1}]

Feladat

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)

(2/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 köbös 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

Gram-Schmidt

P[elem_,l_]:=Module[{},loc=elem-(∑_(i=1)^(Length[l]))(((∫_(-1)^1)elem l[[i]]\[DifferentialD]x)/(∫_(-1)^1)((l[[i]])^2)\[DifferentialD]x)l[[i]];Simplify[(1/loc/.x->1) loc]]

P[x,{1}]

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

P[x^3,{1,x,1/2(3x^2-1)}]

LegendreP[3,x]

P[x^4,{1,x,1/2(3x^2-1),(1/2) (-3 x+5 (x^3))}]

LegendreP[4,x]

L5=P[x^5,{1,x,1/2(3x^2-1),(1/2) (-3 x+5 (x^3)),(1/8) (3-30 (x^2)+35 (x^4))}]

(∫_(-1)^1)L5 1/2(3x^2-1)\[DifferentialD]x==0