Trapezoidal quadrature weights from nodes

11 次查看(过去 30 天)
How can I obtain the trapezoidal quadrature weights from given nodes in MATLAB?

采纳的回答

John D'Errico
John D'Errico 2023-4-10
编辑:John D'Errico 2023-4-11
Hmm. Why not just call trapz? (A suspicious mind here, as this is too easy a question, and there are few good reasons why you need this. But, I can think of at least one. And you already have one answer.)
Simple enough. Given a vector of nodes...
Note that I lack your nodes, so I need to make something up, so I can show how it works.
nodes = [0,sort(rand(1,10)),1];
I chose the very first node at 0 and the last at 1, so I can know the exact result to compare to, and it is easy to compute.
What are the weights for that set of nodes? You do that by composing the weight for each single interval. Then what are the weights for one interval? The area of a trapezoid is simple. You use the average of the function values at the endpoints of the interval, multiplied by the width of the interval. Now put each interval together. It would look like this:
I = (f(1) + f(2))*dx(1)/2 + (f(2) + f(3))*dx(2)/2 + ...
Note that the first and last data points only live in one of the intervals. But you want only the weights, so the coefficients of the function values. This makes the problem pretty simple.
dx = diff(nodes);
W = ([0,dx] + [dx,0])/2;
For this set of nodes, the weights are:
W
W = 1×12
0.0333 0.0626 0.0602 0.0734 0.0586 0.0271 0.0853 0.1982 0.1364 0.0514 0.1262 0.0873
Are those the correct nodes? TRY IT!!!!! For example, I'll pick the function sin(x) on that interval. We know the integral should be 1-cos(1).
fun = @sin;
syms X % verify my claim
int(fun(X),[0,1])
ans = 
Whew! I was right. The exact result (in double precision)
format long g
1 - cos(1)
ans =
0.45969769413186
Now cpompare to trapz.
trapz(nodes,fun(nodes))
ans =
0.458479555032643
Finally, using those nodes is trivially easy now. It is just a dot product.
dot(W,fun(nodes))
ans =
0.458479555032643
As you can see, this replicates the result from trapz, down to the last digit reported.

更多回答(1 个)

Torsten
Torsten 2023-4-10
编辑:Torsten 2023-4-11
Both quadrature weights are 0.5, aren't they ?
1/(x1-x0) * integral_{x=x0}^{x=x1} f(x) dx ~ 0.5*f(x0) + 0.5*f(x1)

类别

Help CenterFile Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by