Trapezoidal quadrature weights from nodes
11 次查看(过去 30 天)
显示 更早的评论
How can I obtain the trapezoidal quadrature weights from given nodes in MATLAB?
0 个评论
采纳的回答
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
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])
Whew! I was right. The exact result (in double precision)
format long g
1 - cos(1)
Now cpompare to trapz.
trapz(nodes,fun(nodes))
Finally, using those nodes is trivially easy now. It is just a dot product.
dot(W,fun(nodes))
As you can see, this replicates the result from trapz, down to the last digit reported.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!