Gauss-Quadratur

05.01.2009, 13:14

tigerbine

Gauss-Quadratur

Kurz gesagt geht man ja hier so vor, dass man die Knoten als Nullstellen von orthogonalen Polynomen bestimmt und dann die zugehörigen Gewichte. Das könnte man dann theoretisch über Integration der Lagrange-Polynome machen. Soweit so gut.

Nun habe ich den leichtsinnigen Weg in ein konkretes Beispiel gewagt... und schwups, ist es da. Das Problem. In meinen Büchern gleichen sich zu berechnenden Integrale und Skalaprodukte dahingehend, dass sie über die gleichen Integrationsgrenzen verfügen...

Theorie

$\begin{eqnarray*} <u,v>:=\int_{a}^{b}\gamma(x)u(x)v(x)~dx \end{eqnarray*}$

Dabei bezeichnet gamma die Gewichtsfunktion. Bzgl. dieses Skalarprodukts orthogonalisert man dann die Monombasis 1,x,x²,... der gewünschten Länge z.B. über Gram-Schmidt oder verwendet das Verfahren von Golub Welsch.

Dabei war das Ziel, folgendes zu berechnen:

$\begin{eqnarray*} \int_a^b~\gamma(x)f(x)~dx \end{eqnarray*}$

Und man erhielt die Gewichte, indem man statt f das IP durch die Knoten integrierte (Lagrange-Darstellung)

$\begin{eqnarray*} \omega_j=\int_a^b\gamma(x)l_j(x)~dx \end{eqnarray*}$

Und die Quadraturformal lautete dann:

$\begin{eqnarray*} I_n(f)=\sum_{j=0}^n~\omega_jf(x_j) \end{eqnarray*}$

Problem:

Nun weichen aber a und b beim Skalarprodukt und dem zu berechnenden Integral ab.

$\begin{eqnarray*} <u,v>=\int_{0}^{1} x^{-\frac{1}{3}}u(x)v(x)~dx \end{eqnarray*}$

Die Gewichtsfunktion lautet also:

$\begin{eqnarray*} \gamma(x)=x^{-\frac{1}{3}} \end{eqnarray*}$

Gesucht ist nun aber das uneigentliche Integral

$\begin{eqnarray*} \int_{1}^{\infty} t^{-\frac{5}{3}}\sin(\frac{1}{t})~dt = \int_{1}^{\infty} t^{-\frac{1}{3}} \cdot t^{-\frac{4}{3}} \cdot \sin(\frac{1}{t})~dt \end{eqnarray*}$ (*)

Damit wäre dann

$\begin{eqnarray*} f(x)=t^{-\frac{4}{3}} \cdot \sin(\frac{1}{t}) \end{eqnarray*}$

Man soll nun (*) mittels Gauss-Quadratur für n=2 (Knoten) approximieren. Die Knoten sind dann die Nullstellen der Polynome, aber wie lauten hier die Gewichte? verwirrt

05.01.2009, 17:08

tigerbine

Auf diesen Beitrag antworten »

RE: Gauss-Quadratur
Also ich habe nun mal Golub Welsh durchgerechnet, das liefert:

$\begin{eqnarray*} q_0(x)=1 \end{eqnarray*}$

$\begin{eqnarray*} q_1(x)=x-\frac{2}{5} \end{eqnarray*}$

$\begin{eqnarray*} q_2(x)=x^2-\frac{10}{11}x+\frac{5}{44} \end{eqnarray*}$

Also die gleichen Polynome wie bei GramSchmidt. Über die EW und EV der Matrix bekomme ich dann die Knoten

$\begin{eqnarray*} x_{1,2}=\frac{5}{11} \pm \frac{3}{22}\sqrt{5} \end{eqnarray*}$

was auch die Nullstellen von q2 sind. Somit sind auch hier beide Verfahren gleich. Die Gewichte ergeben sich gerundet mit GW zu

$\begin{eqnarray*} \omega_1=0.8842 \end{eqnarray*}$

$\begin{eqnarray*} \omega_2=0.6158 \end{eqnarray*}$

In der Summe macht das 3/2, was auch mit der Theorie konform geht (Integration von q0 muss exakt sein).

$\begin{eqnarray*} \int_0^1~x^{-\frac{1}{3}}~dx = \frac{3}{2} \stackrel{!} = \omega_1 + \omega_2 \end{eqnarray*}$

Geht man nun über die Lagrange-Polynome, so ergeben sich, wenn man über [0,1] integriert, die gleichen Werte.

Wie soll man nun aber das uneigentliche Integral approximieren? Wie ist es generell, also wenn es ein bestimmtes Integral wäre?

Würde man auf [0,1] die Knoten/Gewichte bestimmen und dann eine Koordinatentransformation machen?

05.01.2009, 19:37

tigerbine

Auf diesen Beitrag antworten »

RE: Gauss-Quadratur
Ich habe die Transformationsidee mal verfolgt. Da man das uneingentliche Integral ja eigentlich mit Granzwert schreiben müßte, dachte ich man verschiebt die obere Grenze b "langsam" gegen unendlich. Das hat bei mir aber nur die Folge, dass der Wert der Quadraturformel gegen 0 geht.... verwirrt

Würde es nicht auch mehr Sinn machen, wenn man nun an dem Integral interessiert ist,

$\begin{eqnarray*} \int_{1}^{\infty} t^{-\frac{5}{3}}\sin(\frac{1}{t})~dt \end{eqnarray*}$

eine andere Gewichtsfunktion zu nehmen? verwirrt

Oder der Schlüssel liegt in der umgekehrten Transformation... verwirrt

code:

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:

%test 9.1

n=input('b eingeben: ')

disp(' b         Ib        ')
disp('=====================')
for b=2:n
    
    w1=0.8842*(b-1);
    w2=0.6158*(b-1);
    
    x1=1+0.1496*(b-1);
    x2=1+0.7595*(b-1);
    
    fx1=x1^(-4/3)*sin(1/x1);
    fx2=x2^(-4/3)*sin(1/x2);
    
    I=w1*fx1+w2*fx2;
    
    fprintf('  %g      %g   \n',b,I)
end

code:

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:

test91
b eingeben: 100
n =
   100
 b         Ib        
=====================
  2      0.717154   
  3      1.00704   
  4      1.14415   
  5      1.20214   
  6      1.21542   
  7      1.2029   
  8      1.17557   
  9      1.14007   
  10      1.10043   
  11      1.05911   
  12      1.01765   
  13      0.976971   
  14      0.937611   
  15      0.899876   
  16      0.863916   
  17      0.829781   
  18      0.797461   
  19      0.766905   
  20      0.738044   
  21      0.710793   
  22      0.685062   
  23      0.660761   
  24      0.6378   
  25      0.616093   
  26      0.595557   
  27      0.576116   
  28      0.557696   
  29      0.540228   
  30      0.523651   
  31      0.507905   
  32      0.492936   
  33      0.478694   
  34      0.465131   
  35      0.452205   
  36      0.439876   
  37      0.428106   
  38      0.416863   
  39      0.406113   
  40      0.395827   
  41      0.385979   
  42      0.376543   
  43      0.367495   
  44      0.358814   
  45      0.350479   
  46      0.342471   
  47      0.334773   
  48      0.327368   
  49      0.320241   
  50      0.313377   
  51      0.306763   
  52      0.300386   
  53      0.294235   
  54      0.288298   
  55      0.282566   
  56      0.277028   
  57      0.271675   
  58      0.266499   
  59      0.261491   
  60      0.256644   
  61      0.25195   
  62      0.247404   
  63      0.242998   
  64      0.238726   
  65      0.234583   
  66      0.230564   
  67      0.226662   
  68      0.222874   
  69      0.219194   
  70      0.215619   
  71      0.212144   
  72      0.208765   
  73      0.205478   
  74      0.202281   
  75      0.199169   
  76      0.196139   
  77      0.193189   
  78      0.190315   
  79      0.187515   
  80      0.184785   
  81      0.182124   
  82      0.17953   
  83      0.176998   
  84      0.174529   
  85      0.172119   
  86      0.169766   
  87      0.167468   
  88      0.165225   
  89      0.163033   
  90      0.160891   
  91      0.158798   
  92      0.156752   
  93      0.154751   
  94      0.152795   
  95      0.150881   
  96      0.149009   
  97      0.147177   
  98      0.145384   
  99      0.143628   
  100      0.141909

05.01.2009, 23:45

tigerbine

Auf diesen Beitrag antworten »

Ich hab das nun mal mit der anderen Substitution versucht. Kann das mal jemand nachrechnen?

Zitat:

Teilaufgabe a

$\begin{eqnarray*} u_1&&=1 \end{eqnarray*}$

$\begin{eqnarray*} u_2&&= v_2 - \frac{\langle v_2, u_1 \rangle}{\langle u_1,u_1 \rangle} \cdot u_1 =x-\frac{\int_0^1x^{\frac{2}{3}}dx}{\int_0^1x^{-\frac{1}{3}}dx} = x-\frac{\frac{3}{5}}{\frac{3}{2}} \\&&= x-\frac{2}{5} \end{eqnarray*}$

$\begin{eqnarray*} u_3 &&=v_3 - \frac{\langle v_3, u_1 \rangle}{\langle u_1,u_1 \rangle} \cdot u_1 - \frac{\langle v_3, u_2 \rangle}{\langle u_2,u_2 \rangle} \cdot u_2 \\&&=x^2- \frac{10}{11}x+\frac{5}{44} \end{eqnarray*}$

Teilaufgabe b1

Die Knoten der Gaussquadratur sind die Nullstellen von u_3. Da es hier ein quadr. Polynom ist, zieht man die Lösungsformel der Eigenwertbestimmung nach Golub& Welsh vor.

$\begin{eqnarray*} x_{1} =\frac{5}{11} - \frac{3}{22}\sqrt{5} \approx 0.149627 \end{eqnarray*}$

$\begin{eqnarray*} x_{2} =\frac{5}{11} + \frac{3}{22}\sqrt{5} \approx 0.759463 \end{eqnarray*}$

Die Gewichte erhält man über die Gwichtete Integration der Lagrange-Polynome:

$\begin{eqnarray*} \omega_1 = \int_0^1 x^{-\frac{1}{3}}\frac{x-x_2}{x_1-x_2}\dx \approx 0.884164 \end{eqnarray*}$

$\begin{eqnarray*} \omega_2 = \int_0^1 x^{-\frac{1}{3}}\frac{x-x_1}{x_2-x_1}\dx \approx 0.615836 \end{eqnarray*}$

Somit gilt die Approximation:

$\begin{eqnarray*} \int_{0}^{1} x^{-\frac{1}{3}}f(x)~dx \approx w_1f(x_1)+w_2f(x_2) \end{eqnarray*}$

Teilaufgabe b2

Wir interessieren uns nun für das Integral

$\begin{eqnarray*} \int_{1}^{\infty} t^{-\frac{5}{3}}\sin(\frac{1}{t})~dt \end{eqnarray*}$

das auf den ersten Blick nicht zu unserer Vorarbeit passt. Die Verwendung der Variable t gibt uns den entscheidenden Tipp für die Substitution

$\begin{eqnarray*} \int_{\varphi(0)}^{\varphi(1)}~f(t)~dt = \int_{0}^{1}~(f \circ \varphi(x)) \cdot \varphi'(x)~dx \end{eqnarray*}$

$\begin{eqnarray*} \varphi:~[0,1] \to [1,\infty),\quad \varphi(x):=\frac{1}{x}, \frac{1}{t}=x \end{eqnarray*}$

Und wir erhalten schließlich:

$\begin{eqnarray*} \int_{1}^{\infty}~f(t)~dt = \int_{0}^{1}~\frac{1}{x^2}f(\frac{1}{x}) ~dx \end{eqnarray*}$

$\begin{eqnarray*} \int_{1}^{\infty} t^{-\frac{5}{3}}\sin(\frac{1}{t})~dt &&= \int_{1}^{\infty} (\frac{1}{t})^{\frac{5}{3}}\sin(\frac{1}{t})~dt \\ &&=\int_{0}^{1}~\frac{1}{x^2} \cdot x^{\frac{5}{3}}\cdot \sin(x) ~dx \\&&=\int_{0}^{1}~x^{-\frac{1}{3}}\cdot \sin(x) ~dx \end{eqnarray*}$

Somit ergibt sich für die Quadraturformel:

$\begin{eqnarray*} \int_{0}^{1} x^{-\frac{1}{3}}f(x)~dx && \approx w_1f(x_1)+w_2f(x_2) \\&& = w_1 \sin(x_1) + w_2\sin(x_2) \\&& \approx 0.4463 \end{eqnarray*}$

Neue Frage »

Antworten »

Gauss-Quadratur

Verwandte Themen