Mixing
Gaussian quadrature for arbitrary weight function except for the Method of Undetermined Coefficients
The question is as follow:
For the integration
$$int_{-infty}^infty w(x)f(x)dx,$$
where $w(x)$ is any form of distribution, if I want to solve the integration with Gaussian-Hermite quadrature by using look-up table method to get the nodes and weights, how should I do it?
I know the Method of Undetermined Coefficients and Moment-Matching method, but the two methods are troublesome to solve a group of non-linear equation sets. I have ever read that the integration above can be done by the way that the Rackwitz-Fiessler transformation is first used to transform the variable to standard normal variable and then the nodes and weights are derived by using Gaussian-Hermite quadrature. However, I don't know the specific procedures in it.
I really hope somebody can help me out. Thanks!
integration numerical-methods
edited Mar 20 '17 at 14:21
Ian
67.3k25387
asked Mar 20 '17 at 13:59
Lang Wu
85
add a comment |
The question is as follow:
For the integration
$$int_{-infty}^infty w(x)f(x)dx,$$
where $w(x)$ is any form of distribution, if I want to solve the integration with Gaussian-Hermite quadrature by using look-up table method to get the nodes and weights, how should I do it?
I know the Method of Undetermined Coefficients and Moment-Matching method, but the two methods are troublesome to solve a group of non-linear equation sets. I have ever read that the integration above can be done by the way that the Rackwitz-Fiessler transformation is first used to transform the variable to standard normal variable and then the nodes and weights are derived by using Gaussian-Hermite quadrature. However, I don't know the specific procedures in it.
I really hope somebody can help me out. Thanks!
integration numerical-methods
edited Mar 20 '17 at 14:21
Ian
67.3k25387
asked Mar 20 '17 at 13:59
Lang Wu
85
Are you asking "how can I calculate the nodes and weights for a Gaussian quadrature method with an arbitrary nonnegative integrable weight function?" This requires you to compute some integrals of polynomials against $w$ and then to combine these results in some manner. A standard way to do this "combining" step is the Golub-Welsch algorithm, which requires you to solve a tridiagonal eigenproblem, which can usually be done quite quickly for problems of reasonable size.
– Ian
Mar 20 '17 at 14:08
Note that in particular, this tridiagonal eigenproblem in general cannot be solved in closed form (essentially because of the Abel-Ruffini theorem). But this does not really matter; the difficulty is more likely to be in computing the integrals rather than solving the eigenproblem (unless you're trying to make a method with a lot of nodes).
– Ian
Mar 20 '17 at 14:30
Actually, I know one method by searching for orthoganal polynomials but I think it's too confusing. The method you referred to also obscure me. Anyway, thanks for your answer!
– Lang Wu
Mar 20 '17 at 14:32
Alright, I'll give some details in an answer then.
– Ian
Mar 20 '17 at 14:33
I don't need a closed-form solution. You can also choose to introduce me a link or paper about the method.
– Lang Wu
Mar 20 '17 at 14:38
add a comment |
The question is as follow:
For the integration
$$int_{-infty}^infty w(x)f(x)dx,$$
where $w(x)$ is any form of distribution, if I want to solve the integration with Gaussian-Hermite quadrature by using look-up table method to get the nodes and weights, how should I do it?
I know the Method of Undetermined Coefficients and Moment-Matching method, but the two methods are troublesome to solve a group of non-linear equation sets. I have ever read that the integration above can be done by the way that the Rackwitz-Fiessler transformation is first used to transform the variable to standard normal variable and then the nodes and weights are derived by using Gaussian-Hermite quadrature. However, I don't know the specific procedures in it.
I really hope somebody can help me out. Thanks!
integration numerical-methods
edited Mar 20 '17 at 14:21
Ian
67.3k25387
asked Mar 20 '17 at 13:59
Lang Wu
85
The question is as follow:
For the integration
$$int_{-infty}^infty w(x)f(x)dx,$$
where $w(x)$ is any form of distribution, if I want to solve the integration with Gaussian-Hermite quadrature by using look-up table method to get the nodes and weights, how should I do it?
I know the Method of Undetermined Coefficients and Moment-Matching method, but the two methods are troublesome to solve a group of non-linear equation sets. I have ever read that the integration above can be done by the way that the Rackwitz-Fiessler transformation is first used to transform the variable to standard normal variable and then the nodes and weights are derived by using Gaussian-Hermite quadrature. However, I don't know the specific procedures in it.
I really hope somebody can help me out. Thanks!
integration numerical-methods
integration numerical-methods
edited Mar 20 '17 at 14:21
Ian
67.3k25387
asked Mar 20 '17 at 13:59
Lang Wu
85
edited Mar 20 '17 at 14:21
Ian
67.3k25387
asked Mar 20 '17 at 13:59
Lang Wu
85
edited Mar 20 '17 at 14:21
Ian
67.3k25387
edited Mar 20 '17 at 14:21
Ian
67.3k25387
edited Mar 20 '17 at 14:21
Ian
67.3k25387
67.3k25387
asked Mar 20 '17 at 13:59
Lang Wu
85
asked Mar 20 '17 at 13:59
Lang Wu
85
asked Mar 20 '17 at 13:59
Lang Wu
85
85
Are you asking "how can I calculate the nodes and weights for a Gaussian quadrature method with an arbitrary nonnegative integrable weight function?" This requires you to compute some integrals of polynomials against $w$ and then to combine these results in some manner. A standard way to do this "combining" step is the Golub-Welsch algorithm, which requires you to solve a tridiagonal eigenproblem, which can usually be done quite quickly for problems of reasonable size.
– Ian
Mar 20 '17 at 14:08
Note that in particular, this tridiagonal eigenproblem in general cannot be solved in closed form (essentially because of the Abel-Ruffini theorem). But this does not really matter; the difficulty is more likely to be in computing the integrals rather than solving the eigenproblem (unless you're trying to make a method with a lot of nodes).
– Ian
Mar 20 '17 at 14:30
Actually, I know one method by searching for orthoganal polynomials but I think it's too confusing. The method you referred to also obscure me. Anyway, thanks for your answer!
– Lang Wu
Mar 20 '17 at 14:32
Alright, I'll give some details in an answer then.
– Ian
Mar 20 '17 at 14:33
I don't need a closed-form solution. You can also choose to introduce me a link or paper about the method.
– Lang Wu
Mar 20 '17 at 14:38
add a comment |
Are you asking "how can I calculate the nodes and weights for a Gaussian quadrature method with an arbitrary nonnegative integrable weight function?" This requires you to compute some integrals of polynomials against $w$ and then to combine these results in some manner. A standard way to do this "combining" step is the Golub-Welsch algorithm, which requires you to solve a tridiagonal eigenproblem, which can usually be done quite quickly for problems of reasonable size.
– Ian
Mar 20 '17 at 14:08
Note that in particular, this tridiagonal eigenproblem in general cannot be solved in closed form (essentially because of the Abel-Ruffini theorem). But this does not really matter; the difficulty is more likely to be in computing the integrals rather than solving the eigenproblem (unless you're trying to make a method with a lot of nodes).
– Ian
Mar 20 '17 at 14:30
Actually, I know one method by searching for orthoganal polynomials but I think it's too confusing. The method you referred to also obscure me. Anyway, thanks for your answer!
– Lang Wu
Mar 20 '17 at 14:32
Alright, I'll give some details in an answer then.
– Ian
Mar 20 '17 at 14:33
I don't need a closed-form solution. You can also choose to introduce me a link or paper about the method.
– Lang Wu
Mar 20 '17 at 14:38
Are you asking "how can I calculate the nodes and weights for a Gaussian quadrature method with an arbitrary nonnegative integrable weight function?" This requires you to compute some integrals of polynomials against $w$ and then to combine these results in some manner. A standard way to do this "combining" step is the Golub-Welsch algorithm, which requires you to solve a tridiagonal eigenproblem, which can usually be done quite quickly for problems of reasonable size.
– Ian
Mar 20 '17 at 14:08
Are you asking "how can I calculate the nodes and weights for a Gaussian quadrature method with an arbitrary nonnegative integrable weight function?" This requires you to compute some integrals of polynomials against $w$ and then to combine these results in some manner. A standard way to do this "combining" step is the Golub-Welsch algorithm, which requires you to solve a tridiagonal eigenproblem, which can usually be done quite quickly for problems of reasonable size.
– Ian
Mar 20 '17 at 14:08
Note that in particular, this tridiagonal eigenproblem in general cannot be solved in closed form (essentially because of the Abel-Ruffini theorem). But this does not really matter; the difficulty is more likely to be in computing the integrals rather than solving the eigenproblem (unless you're trying to make a method with a lot of nodes).
– Ian
Mar 20 '17 at 14:30
Note that in particular, this tridiagonal eigenproblem in general cannot be solved in closed form (essentially because of the Abel-Ruffini theorem). But this does not really matter; the difficulty is more likely to be in computing the integrals rather than solving the eigenproblem (unless you're trying to make a method with a lot of nodes).
– Ian
Mar 20 '17 at 14:30
Actually, I know one method by searching for orthoganal polynomials but I think it's too confusing. The method you referred to also obscure me. Anyway, thanks for your answer!
– Lang Wu
Mar 20 '17 at 14:32
Actually, I know one method by searching for orthoganal polynomials but I think it's too confusing. The method you referred to also obscure me. Anyway, thanks for your answer!
– Lang Wu
Mar 20 '17 at 14:32
Alright, I'll give some details in an answer then.
– Ian
Mar 20 '17 at 14:33
Alright, I'll give some details in an answer then.
– Ian
Mar 20 '17 at 14:33
I don't need a closed-form solution. You can also choose to introduce me a link or paper about the method.
– Lang Wu
Mar 20 '17 at 14:38
I don't need a closed-form solution. You can also choose to introduce me a link or paper about the method.
– Lang Wu
Mar 20 '17 at 14:38
add a comment |
1 Answer
1
active
oldest
votes
The abbreviated explanation for implementation of Gaussian quadrature looks like this. I'll define:
- the inner product $(f,g)=int_{-infty}^infty w(x) f(x) g(x) dx$ on suitable functions $f,g in L^2_w$ (if you don't know what $L^2_w$ means, ignore it).
- the corresponding norm $| f |=(f,f)^{1/2}$.
Next I'll define two sequences of polynomials $p_k,q_k$ with
- $p_k=q_k/| q_k |$
$q_0=1$- $q_1=x-(xp_0,p_0)p_0$
- $q_k=xp_{k-1}-(xp_{k-1},p_{k-1})p_{k-1}-(xp_{k-1},p_{k-2})p_{k-2},k geq 2$
This last equation is called a "three term recurrence relation". In this setting, $p_k$ has degree $k$, has norm $1$, and is orthogonal to all the other $p_j$.
Define $a_k=(xp_k,p_k)$ for $k geq 0$, and $b_k=(xp_k,p_{k-1})$ for $k geq 1$. (Note that these are exactly the integrals computed to generate the $p_k$.) Then we define the $n times n$ matrix $T_{ij}$ with $T_{ii}=a_{i-1}$ for $i geq 1$ and $T_{i-1,i}=T_{i,i-1}=(b_{i-1})^{1/2}$ for $i geq 2$ and all other $T_{ij}=0$. Finally we define $c=(1,1)$ or in other words $c=int_{-infty}^infty w(x) dx$. ($c$ had to come in somewhere, because the $p_k$ remain invariant under multiplication of $w$ by a constant.)
Then the nodes $x_i$ of the $n$ point Gaussian quadrature method are the eigenvalues of $T$. The weight $w_i = c (y^{(i)}_1)^2$ where $y^{(i)}$ is a unit eigenvector with eigenvalue $x_i$. Since this is a symmetric tridiagonal eigenproblem, it is pretty easy to solve, relative to its size. The task of computing the $a_i,b_i$ will probably take more time than solving the eigenproblem (unless they can be computed in closed form, which they can for certain simple $w$).
This procedure is called the Golub-Welsch algorithm. It works primarily because of the fundamental theorem of Gaussian quadrature, which says that the nodes of a $n$ point Gaussian quadrature are the roots of $p_n$. You can find more details at https://en.wikipedia.org/wiki/Gaussian_quadrature
I'll add that technically this is equivalent to "undetermined coefficients", i.e. to setting up the nonlinear system $sum_{i=1}^n w_i x_i^k = int_{-infty}^infty w(x) x^k dx$ for $k=0,1,dots,2n-1$. But in practice Golub-Welsch will be much easier than solving that nonlinear system. That is because:
- the nonlinear system is poorly scaled. (It would be well-scaled if you used $p_k$ instead of $x^k$.)
- Casting the nonlinear system into a symmetric tridiagonal eigenproblem makes it far easier.
answered Mar 20 '17 at 15:18
Ian
67.3k25387
It would be helpful. Thank you very much!
– Lang Wu
Mar 21 '17 at 13:57
@LangWu What would?
– Ian
Mar 21 '17 at 14:02
I mean I need time to understand this method and try to use it since I haven't systematically learned the Gaussian quadrature but just read it from Wikipedia. On the other hand, this question is just a part of what I want to solve. Maybe I will put forward it before long and sincerely hope you can help me then.
– Lang Wu
Mar 21 '17 at 14:24
add a comment |
Your Answer
StackExchange.ifUsing("editor", function () {
return StackExchange.using("mathjaxEditing", function () {
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
});
});
}, "mathjax-editing");
StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "69"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);
StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});
function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: true,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
noCode: true, onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});
}
});
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f2195054%2fgaussian-quadrature-for-arbitrary-weight-function-except-for-the-method-of-undet%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
1 Answer
1
active
oldest
votes
1 Answer
1
active
oldest
votes
active
oldest
votes
active
oldest
votes
The abbreviated explanation for implementation of Gaussian quadrature looks like this. I'll define:
- the inner product $(f,g)=int_{-infty}^infty w(x) f(x) g(x) dx$ on suitable functions $f,g in L^2_w$ (if you don't know what $L^2_w$ means, ignore it).
- the corresponding norm $| f |=(f,f)^{1/2}$.
Next I'll define two sequences of polynomials $p_k,q_k$ with
- $p_k=q_k/| q_k |$
$q_0=1$- $q_1=x-(xp_0,p_0)p_0$
- $q_k=xp_{k-1}-(xp_{k-1},p_{k-1})p_{k-1}-(xp_{k-1},p_{k-2})p_{k-2},k geq 2$
This last equation is called a "three term recurrence relation". In this setting, $p_k$ has degree $k$, has norm $1$, and is orthogonal to all the other $p_j$.
Define $a_k=(xp_k,p_k)$ for $k geq 0$, and $b_k=(xp_k,p_{k-1})$ for $k geq 1$. (Note that these are exactly the integrals computed to generate the $p_k$.) Then we define the $n times n$ matrix $T_{ij}$ with $T_{ii}=a_{i-1}$ for $i geq 1$ and $T_{i-1,i}=T_{i,i-1}=(b_{i-1})^{1/2}$ for $i geq 2$ and all other $T_{ij}=0$. Finally we define $c=(1,1)$ or in other words $c=int_{-infty}^infty w(x) dx$. ($c$ had to come in somewhere, because the $p_k$ remain invariant under multiplication of $w$ by a constant.)
Then the nodes $x_i$ of the $n$ point Gaussian quadrature method are the eigenvalues of $T$. The weight $w_i = c (y^{(i)}_1)^2$ where $y^{(i)}$ is a unit eigenvector with eigenvalue $x_i$. Since this is a symmetric tridiagonal eigenproblem, it is pretty easy to solve, relative to its size. The task of computing the $a_i,b_i$ will probably take more time than solving the eigenproblem (unless they can be computed in closed form, which they can for certain simple $w$).
This procedure is called the Golub-Welsch algorithm. It works primarily because of the fundamental theorem of Gaussian quadrature, which says that the nodes of a $n$ point Gaussian quadrature are the roots of $p_n$. You can find more details at https://en.wikipedia.org/wiki/Gaussian_quadrature
I'll add that technically this is equivalent to "undetermined coefficients", i.e. to setting up the nonlinear system $sum_{i=1}^n w_i x_i^k = int_{-infty}^infty w(x) x^k dx$ for $k=0,1,dots,2n-1$. But in practice Golub-Welsch will be much easier than solving that nonlinear system. That is because:
- the nonlinear system is poorly scaled. (It would be well-scaled if you used $p_k$ instead of $x^k$.)
- Casting the nonlinear system into a symmetric tridiagonal eigenproblem makes it far easier.
answered Mar 20 '17 at 15:18
Ian
67.3k25387
It would be helpful. Thank you very much!
– Lang Wu
Mar 21 '17 at 13:57
@LangWu What would?
– Ian
Mar 21 '17 at 14:02
I mean I need time to understand this method and try to use it since I haven't systematically learned the Gaussian quadrature but just read it from Wikipedia. On the other hand, this question is just a part of what I want to solve. Maybe I will put forward it before long and sincerely hope you can help me then.
– Lang Wu
Mar 21 '17 at 14:24
add a comment |
The abbreviated explanation for implementation of Gaussian quadrature looks like this. I'll define:
- the inner product $(f,g)=int_{-infty}^infty w(x) f(x) g(x) dx$ on suitable functions $f,g in L^2_w$ (if you don't know what $L^2_w$ means, ignore it).
- the corresponding norm $| f |=(f,f)^{1/2}$.
Next I'll define two sequences of polynomials $p_k,q_k$ with
- $p_k=q_k/| q_k |$
$q_0=1$- $q_1=x-(xp_0,p_0)p_0$
- $q_k=xp_{k-1}-(xp_{k-1},p_{k-1})p_{k-1}-(xp_{k-1},p_{k-2})p_{k-2},k geq 2$
This last equation is called a "three term recurrence relation". In this setting, $p_k$ has degree $k$, has norm $1$, and is orthogonal to all the other $p_j$.
Define $a_k=(xp_k,p_k)$ for $k geq 0$, and $b_k=(xp_k,p_{k-1})$ for $k geq 1$. (Note that these are exactly the integrals computed to generate the $p_k$.) Then we define the $n times n$ matrix $T_{ij}$ with $T_{ii}=a_{i-1}$ for $i geq 1$ and $T_{i-1,i}=T_{i,i-1}=(b_{i-1})^{1/2}$ for $i geq 2$ and all other $T_{ij}=0$. Finally we define $c=(1,1)$ or in other words $c=int_{-infty}^infty w(x) dx$. ($c$ had to come in somewhere, because the $p_k$ remain invariant under multiplication of $w$ by a constant.)
Then the nodes $x_i$ of the $n$ point Gaussian quadrature method are the eigenvalues of $T$. The weight $w_i = c (y^{(i)}_1)^2$ where $y^{(i)}$ is a unit eigenvector with eigenvalue $x_i$. Since this is a symmetric tridiagonal eigenproblem, it is pretty easy to solve, relative to its size. The task of computing the $a_i,b_i$ will probably take more time than solving the eigenproblem (unless they can be computed in closed form, which they can for certain simple $w$).
This procedure is called the Golub-Welsch algorithm. It works primarily because of the fundamental theorem of Gaussian quadrature, which says that the nodes of a $n$ point Gaussian quadrature are the roots of $p_n$. You can find more details at https://en.wikipedia.org/wiki/Gaussian_quadrature
I'll add that technically this is equivalent to "undetermined coefficients", i.e. to setting up the nonlinear system $sum_{i=1}^n w_i x_i^k = int_{-infty}^infty w(x) x^k dx$ for $k=0,1,dots,2n-1$. But in practice Golub-Welsch will be much easier than solving that nonlinear system. That is because:
- the nonlinear system is poorly scaled. (It would be well-scaled if you used $p_k$ instead of $x^k$.)
- Casting the nonlinear system into a symmetric tridiagonal eigenproblem makes it far easier.
answered Mar 20 '17 at 15:18
Ian
67.3k25387
It would be helpful. Thank you very much!
– Lang Wu
Mar 21 '17 at 13:57
@LangWu What would?
– Ian
Mar 21 '17 at 14:02
I mean I need time to understand this method and try to use it since I haven't systematically learned the Gaussian quadrature but just read it from Wikipedia. On the other hand, this question is just a part of what I want to solve. Maybe I will put forward it before long and sincerely hope you can help me then.
– Lang Wu
Mar 21 '17 at 14:24
add a comment |
The abbreviated explanation for implementation of Gaussian quadrature looks like this. I'll define:
- the inner product $(f,g)=int_{-infty}^infty w(x) f(x) g(x) dx$ on suitable functions $f,g in L^2_w$ (if you don't know what $L^2_w$ means, ignore it).
- the corresponding norm $| f |=(f,f)^{1/2}$.
Next I'll define two sequences of polynomials $p_k,q_k$ with
- $p_k=q_k/| q_k |$
$q_0=1$- $q_1=x-(xp_0,p_0)p_0$
- $q_k=xp_{k-1}-(xp_{k-1},p_{k-1})p_{k-1}-(xp_{k-1},p_{k-2})p_{k-2},k geq 2$
This last equation is called a "three term recurrence relation". In this setting, $p_k$ has degree $k$, has norm $1$, and is orthogonal to all the other $p_j$.
Define $a_k=(xp_k,p_k)$ for $k geq 0$, and $b_k=(xp_k,p_{k-1})$ for $k geq 1$. (Note that these are exactly the integrals computed to generate the $p_k$.) Then we define the $n times n$ matrix $T_{ij}$ with $T_{ii}=a_{i-1}$ for $i geq 1$ and $T_{i-1,i}=T_{i,i-1}=(b_{i-1})^{1/2}$ for $i geq 2$ and all other $T_{ij}=0$. Finally we define $c=(1,1)$ or in other words $c=int_{-infty}^infty w(x) dx$. ($c$ had to come in somewhere, because the $p_k$ remain invariant under multiplication of $w$ by a constant.)
Then the nodes $x_i$ of the $n$ point Gaussian quadrature method are the eigenvalues of $T$. The weight $w_i = c (y^{(i)}_1)^2$ where $y^{(i)}$ is a unit eigenvector with eigenvalue $x_i$. Since this is a symmetric tridiagonal eigenproblem, it is pretty easy to solve, relative to its size. The task of computing the $a_i,b_i$ will probably take more time than solving the eigenproblem (unless they can be computed in closed form, which they can for certain simple $w$).
This procedure is called the Golub-Welsch algorithm. It works primarily because of the fundamental theorem of Gaussian quadrature, which says that the nodes of a $n$ point Gaussian quadrature are the roots of $p_n$. You can find more details at https://en.wikipedia.org/wiki/Gaussian_quadrature
I'll add that technically this is equivalent to "undetermined coefficients", i.e. to setting up the nonlinear system $sum_{i=1}^n w_i x_i^k = int_{-infty}^infty w(x) x^k dx$ for $k=0,1,dots,2n-1$. But in practice Golub-Welsch will be much easier than solving that nonlinear system. That is because:
- the nonlinear system is poorly scaled. (It would be well-scaled if you used $p_k$ instead of $x^k$.)
- Casting the nonlinear system into a symmetric tridiagonal eigenproblem makes it far easier.
answered Mar 20 '17 at 15:18
Ian
67.3k25387
The abbreviated explanation for implementation of Gaussian quadrature looks like this. I'll define:
- the inner product $(f,g)=int_{-infty}^infty w(x) f(x) g(x) dx$ on suitable functions $f,g in L^2_w$ (if you don't know what $L^2_w$ means, ignore it).
- the corresponding norm $| f |=(f,f)^{1/2}$.
Next I'll define two sequences of polynomials $p_k,q_k$ with
- $p_k=q_k/| q_k |$
$q_0=1$- $q_1=x-(xp_0,p_0)p_0$
- $q_k=xp_{k-1}-(xp_{k-1},p_{k-1})p_{k-1}-(xp_{k-1},p_{k-2})p_{k-2},k geq 2$
This last equation is called a "three term recurrence relation". In this setting, $p_k$ has degree $k$, has norm $1$, and is orthogonal to all the other $p_j$.
Define $a_k=(xp_k,p_k)$ for $k geq 0$, and $b_k=(xp_k,p_{k-1})$ for $k geq 1$. (Note that these are exactly the integrals computed to generate the $p_k$.) Then we define the $n times n$ matrix $T_{ij}$ with $T_{ii}=a_{i-1}$ for $i geq 1$ and $T_{i-1,i}=T_{i,i-1}=(b_{i-1})^{1/2}$ for $i geq 2$ and all other $T_{ij}=0$. Finally we define $c=(1,1)$ or in other words $c=int_{-infty}^infty w(x) dx$. ($c$ had to come in somewhere, because the $p_k$ remain invariant under multiplication of $w$ by a constant.)
Then the nodes $x_i$ of the $n$ point Gaussian quadrature method are the eigenvalues of $T$. The weight $w_i = c (y^{(i)}_1)^2$ where $y^{(i)}$ is a unit eigenvector with eigenvalue $x_i$. Since this is a symmetric tridiagonal eigenproblem, it is pretty easy to solve, relative to its size. The task of computing the $a_i,b_i$ will probably take more time than solving the eigenproblem (unless they can be computed in closed form, which they can for certain simple $w$).
This procedure is called the Golub-Welsch algorithm. It works primarily because of the fundamental theorem of Gaussian quadrature, which says that the nodes of a $n$ point Gaussian quadrature are the roots of $p_n$. You can find more details at https://en.wikipedia.org/wiki/Gaussian_quadrature
I'll add that technically this is equivalent to "undetermined coefficients", i.e. to setting up the nonlinear system $sum_{i=1}^n w_i x_i^k = int_{-infty}^infty w(x) x^k dx$ for $k=0,1,dots,2n-1$. But in practice Golub-Welsch will be much easier than solving that nonlinear system. That is because:
- the nonlinear system is poorly scaled. (It would be well-scaled if you used $p_k$ instead of $x^k$.)
- Casting the nonlinear system into a symmetric tridiagonal eigenproblem makes it far easier.
answered Mar 20 '17 at 15:18
Ian
67.3k25387
edited Nov 20 '18 at 15:19
answered Mar 20 '17 at 15:18
Ian
67.3k25387
answered Mar 20 '17 at 15:18
Ian
67.3k25387
answered Mar 20 '17 at 15:18
Ian
67.3k25387
67.3k25387
It would be helpful. Thank you very much!
– Lang Wu
Mar 21 '17 at 13:57
@LangWu What would?
– Ian
Mar 21 '17 at 14:02
I mean I need time to understand this method and try to use it since I haven't systematically learned the Gaussian quadrature but just read it from Wikipedia. On the other hand, this question is just a part of what I want to solve. Maybe I will put forward it before long and sincerely hope you can help me then.
– Lang Wu
Mar 21 '17 at 14:24
add a comment |
It would be helpful. Thank you very much!
– Lang Wu
Mar 21 '17 at 13:57
@LangWu What would?
– Ian
Mar 21 '17 at 14:02
I mean I need time to understand this method and try to use it since I haven't systematically learned the Gaussian quadrature but just read it from Wikipedia. On the other hand, this question is just a part of what I want to solve. Maybe I will put forward it before long and sincerely hope you can help me then.
– Lang Wu
Mar 21 '17 at 14:24
It would be helpful. Thank you very much!
– Lang Wu
Mar 21 '17 at 13:57
It would be helpful. Thank you very much!
– Lang Wu
Mar 21 '17 at 13:57
@LangWu What would?
– Ian
Mar 21 '17 at 14:02
@LangWu What would?
– Ian
Mar 21 '17 at 14:02
I mean I need time to understand this method and try to use it since I haven't systematically learned the Gaussian quadrature but just read it from Wikipedia. On the other hand, this question is just a part of what I want to solve. Maybe I will put forward it before long and sincerely hope you can help me then.
– Lang Wu
Mar 21 '17 at 14:24
I mean I need time to understand this method and try to use it since I haven't systematically learned the Gaussian quadrature but just read it from Wikipedia. On the other hand, this question is just a part of what I want to solve. Maybe I will put forward it before long and sincerely hope you can help me then.
– Lang Wu
Mar 21 '17 at 14:24
add a comment |
Thanks for contributing an answer to Mathematics Stack Exchange!
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
Use MathJax to format equations. MathJax reference.
To learn more, see our tips on writing great answers.
Some of your past answers have not been well-received, and you're in danger of being blocked from answering.
Please pay close attention to the following guidance:
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
To learn more, see our tips on writing great answers.
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f2195054%2fgaussian-quadrature-for-arbitrary-weight-function-except-for-the-method-of-undet%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Are you asking "how can I calculate the nodes and weights for a Gaussian quadrature method with an arbitrary nonnegative integrable weight function?" This requires you to compute some integrals of polynomials against $w$ and then to combine these results in some manner. A standard way to do this "combining" step is the Golub-Welsch algorithm, which requires you to solve a tridiagonal eigenproblem, which can usually be done quite quickly for problems of reasonable size.
– Ian
Mar 20 '17 at 14:08
Note that in particular, this tridiagonal eigenproblem in general cannot be solved in closed form (essentially because of the Abel-Ruffini theorem). But this does not really matter; the difficulty is more likely to be in computing the integrals rather than solving the eigenproblem (unless you're trying to make a method with a lot of nodes).
– Ian
Mar 20 '17 at 14:30
Actually, I know one method by searching for orthoganal polynomials but I think it's too confusing. The method you referred to also obscure me. Anyway, thanks for your answer!
– Lang Wu
Mar 20 '17 at 14:32
Alright, I'll give some details in an answer then.
– Ian
Mar 20 '17 at 14:33
I don't need a closed-form solution. You can also choose to introduce me a link or paper about the method.
– Lang Wu
Mar 20 '17 at 14:38