is there a matlab code for gaussion integration!!!
37 views (last 30 days)
Show older comments
Aalaa Abu alrob
on 26 Oct 2018
Edited: Aalaa Abu alrob
on 12 Nov 2018
syms x;
syms E;
syms r;
syms D;
syms c;
syms m;
g(x)=exp((-2./h).*(sqrt(2*m*(E-c))).*(sqrt(x*(x+r))-r*log((sqrt(r)+sqrt(x+r))-D)))
0 Comments
Accepted Answer
John D'Errico
on 26 Oct 2018
Edited: John D'Errico
on 26 Oct 2018
Did you look on the file exchange? (Clearly not.)
They both look decent, though guassquad is purely a gauss-legendre code, gaussg a more general code for standard weight functions, though it is a guass-kronrod scheme.
In fact, even my sympoly toolbox provides gaussian quadrature rules, thus the weights and nodes for any standard weight function class. So look at the guassquadrule function.
So using guassquadrule from sympoly, the nodes and weights for a 4 point Gauss-Hermite rule would be:
[nodes,weights] = gaussquadrule(4,'hermite')
nodes =
-1.65068012388579 -0.52464762327529 0.52464762327529 1.65068012388579
weights =
0.081312835447245 0.804914090005513 0.804914090005513 0.0813128354472448
To compute the integral of (x^5 - 3*x^2 - 2)*exp(-x.^2), from -inf to inf, I might do:
fun = @(x) (x.^5 -3*x.^2 - 2);
dot(weights,fun(nodes))
ans =
-6.2035884781693
Syms confirms that it was correct:
syms x
int((x.^5 -3*x.^2 - 2)*exp(-x^2),[-inf,inf])
ans =
-(7*pi^(1/2))/2
vpa(ans)
ans =
-6.203588478169306095543586191694
However, will you find a tool that performs gaussian integration for a function that contains symbolic parameters? No, you won't find anything for a good reason. Gaussian integration is a numerical integration procedure, not really designed to do symbolic integration.
Can you use gaussian quadrature for a symbolic function? Well, yes.
So, here a Gauss-Laguerre integral of (a*x*2 + b*x + c)*exp(-x), over the domain [0,inf].
syms a b c
fun = @(x) a*x.^2 + b*x + c;
[nodes,weights] = gaussquadrule(5,'laguerre');
vpa(dot(weights,fun(nodes)),13)
ans =
2.0*a + 1.0*b + 1.0*c
What Gaussian quadrature rule you expect to apply to this function:
g(x)=exp((-2./h).*(sqrt(2*m*(E-c))).*(sqrt(x*(x+r))-r*log((sqrt(r)+sqrt(x+r))-D)))
I have no idea, since I don't see any domain provided, nor do I see any standard weight function in there that can be extracted. So I have a funny feeling that you may be confused as to the purpose and utility of classical Gaussian quadratures.
If you look at Gaussian quadrature rules, they presume a weight function from among several standard forms, AND a domain of integration. Gauss-Legendre assumes a unit weight function, so is applicable to integration of a general function, over the interval [-1,1]. For example, suppose you wanted to compute the integral of cos(exp(x)), over the interval [0,1].
syms x
int(cos(exp(x)),[-1,1])
ans =
cosint(exp(1)) - cosint(exp(-1))
vpa(ans)
ans =
0.67038594208938451613294763492665
Luckily, the symbolic TB is smart enough to do the work, because I'm way too sleepy to think.
Gauss-Legendre presumes an interval of [-1,1], although the nodes and weights can be transformed for other intervals.
fun = @(x) cos(exp(x));
[nodes,weights] = gaussquadrule(15,'legendre');
dot(weights,fun(nodes))
ans =
0.670385942089469
6 Comments
John D'Errico
on 1 Nov 2018
CAN you do it then? There are still issues you need to consider.
pretty(vpa(g,5))
34 34
exp(2.6418 10 log(sqrt(x + 2.0) + 1.6923) - 1.3209 10 sqrt(x (x + 2.0)))
We can rewrite this as
g(x) = exp(1.3209e34*H(x))
where H(x) is:
H = @(x) 2*log(sqrt(x+2.0)+1.6923) - sqrt(x.*(x+2.0));
ezplot(H,[-1,1])

Now, you need to recognize that nothing was plotted for negative x. Why not? Because the result is complex for negative x.
H(-.1)
ans =
2.2438 - 0.43589i
Just as bad is the fact that this function is of the general form exp(1.3209e34*H(x)), where H(x) is on the general order of 1.
Do you have any clue how large exp(1.3209e34) is?
By way of comparison, the mass of the sun, measured in grams, is a smaller number.
So, can you compute that integral? Sigh. Does it mean anything?
More Answers (1)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!