Einführung in Mathematica XXIV

Integration mit orthogonalen Polynomen

© Jürgen Schnack, 2012/05/24

Problem

Funktion f[x] soll numerisch integriert werden.
Erzeuge dazu zuerst orthogonales Funktionensystem auf [-1,1].

f[x_]=Exp[Sin[5*x]];

Plot[f[x],{x,-2,2},PlotRange->{-3,3},PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick},{Magenta,Thick},{Cyan,Thick}},BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

NIntegrate[f[x],{x,-1,1}]

Ausgangspunkt Monome (x^n)

NMax=4;

b[x_]=Table[(x^n),{n,0,NMax}]

Plot[{b[x][[1]],b[x][[2]],b[x][[3]],b[x][[4]],b[x][[5]]}, {x, -1, 1},
PlotRange -> {-1, 1},
PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick},{Magenta,Thick},{Cyan,Thick}},BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

Erzeugung der orthogonalen Basis

1. Zustand

Normierung=(∫_(-1)^1)(b[x][[1]]*b[x][[1]])\[DifferentialD]x

onb1[x_]=b[x][[1]]/(Sqrt[Normierung])

Erzeugung der orthogonalen Basis

2. Zustand

Ueberlapp=(∫_(-1)^1)(onb1[x]*b[x][[2]])\[DifferentialD]x

ounb2[x_]=b[x][[2]]-Ueberlapp*onb1[x]

Normierung=(∫_(-1)^1)(ounb2[x]*ounb2[x])\[DifferentialD]x

onb2[x_]=ounb2[x]/(Sqrt[Normierung])

onb[x_]={onb1[x],onb2[x]}

onb[x][[2]]

Erzeugung der orthogonalen Basis

3. Zustand

NZst=3;

Ueber=Table[(∫_(-1)^1)(onb[x][[n]]*b[x][[NZst]])\[DifferentialD]x,{n,1,NZst-1}]

ounbNEU[x_]=b[x][[NZst]]-(∑_(k=1)^(NZst-1))(Ueber[[k]]*onb[x][[k]])

Normierung=(∫_(-1)^1)(ounbNEU[x]*ounbNEU[x])\[DifferentialD]x

onbNEU[x_]=ounbNEU[x]/(Sqrt[Normierung])

onb[x_]=Append[onb[x],onbNEU[x]]

onb[x][[3]]

Darstellung der ersten orthogonalen Polynome

Plot[{onb[x][[1]],onb[x][[2]],onb[x][[3]]},{x,-1,1},PlotRange->{-1,1},PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick},{Magenta,Thick},{Cyan,Thick}},BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

Erzeugung der orthogonalen Basis

alle weiteren Zustände

For[NZst=4,NZst<NMax+2,NZst++,
Ueber = Table[(∫_(-1)^1)(onb[x][[n]]*b[x][[NZst]])\[DifferentialD]x, {n, 1, NZst - 1}];ounbNEU[x_] = b[x][[NZst]] - (∑_(k=1)^(NZst-1))(Ueber[[k]]*onb[x][[k]]);
Normierung = (∫_(-1)^1)(ounbNEU[x]*ounbNEU[x])\[DifferentialD]x;

onbNEU[x_] = ounbNEU[x]/Sqrt[Normierung];
onb[x_] = Append[onb[x], onbNEU[x]]];

onb[x]

Darstellung der orthogonalen Basis

alle Zustände

Plot[{onb[x][[1]],onb[x][[2]],onb[x][[3]],onb[x][[4]],onb[x][[5]]},
{x,-1,1},PlotRange->{-1,1},PlotStyle->{{Red,Thick},{Green,Thick},{Blue,Thick},{Magenta,Thick},{Cyan,Thick}},BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

Nullstellen des Polynoms zur Ordnung NMax

alle weiteren Zustände

Nullstellen=Solve[onb[x][[NMax+1]]==0,x]

TNullstellen=Table[x/.Nullstellen[[k]],{k,1,NMax}]

Berechnung der Gewichte

Lagrange - Polynome

LPolynome[y_]=Table[(∏_(k=1)^NMax)(If[k==n,1,(((y-TNullstellen[[k]]))/((TNullstellen[[n]]-TNullstellen[[k]])))]),{n,1,NMax}]

Berechnung der Gewichte

Gewichte

Wi=Table[(∫_(-1)^1)(LPolynome[y][[n]])\[DifferentialD]y,{n,1,NMax}]

Integral

IntApprox=N[(∑_(k=1)^NMax)(Wi[[k]]*f[TNullstellen[[k]]])]

Nutzung von 32 Nullstellen und Gewichten

Gewichte

G={{-0.997263861849, 0.00701814576495},{-0.985611511545, 0.0162774265831},{-0.964762255588, 0.0253910098329},{-0.934906075938, 0.0342745478477},
{-0.896321155766, 0.0428359896785},{-0.849367613733, 0.0509978738117},{-0.794483795968, 0.0586839394615},{-0.73218211874 ,0.0658220603578},{-0.66304426693 ,0.0723456094297}
,{-0.587715757241, 0.078193695762},{-0.506899908932 ,0.083311711103}
,{-0.421351276131 ,0.0876518688047},{-0.331868602282 ,0.0911736454878}
,{-0.239287362252 ,0.0938441590423},{-0.144471961583 ,0.0956384754512}
,{-0.0483076656877,0.0965398415811} ,{0.0483076656877,0.0965398415811}
 ,{0.144471961583,0.0956384754512} ,{0.239287362252,0.0938441590423}
 ,{0.331868602282,0.0911736454878} ,{0.421351276131,0.0876518688047}
 ,{0.506899908932,0.083311711103} ,{0.587715757241,0.078193695762}
 ,{0.66304426693,0.0723456094297} ,{0.73218211874,0.0658220603578}
 ,{0.794483795968,0.0586839394615} ,{0.849367613733,0.0509978738117}
 ,{0.896321155766,0.0428359896785} ,{0.934906075938,0.0342745478477}
 ,{0.964762255588,0.0253910098329} ,{0.985611511545,0.0162774265831}
 ,{0.997263861849,0.00701814576495}}

Integral

IntApprox=N[(∑_(k=1)^(Length[G]))(G[[k,2]]*f[G[[k,1]]])]

Zum Vergleich: Trapet mit 32+1 Stützstellen

IntTrapez=(Sum[f[-1.0+k*2.0/Length[G]],{k,0,Length[G]}]-0.5*f[-1.0]-0.5*f[1.0])*2.0/Length[G]