forked from rrhiemstra/quadrules
-
Notifications
You must be signed in to change notification settings - Fork 1
Expand file tree
/
Copy pathHowTo.m
More file actions
31 lines (26 loc) · 847 Bytes
/
HowTo.m
File metadata and controls
31 lines (26 loc) · 847 Bytes
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
% initialize
clear; close all; clc;
% choose parameters
p = 3; % choose polynomial degree
r = 0; % regularity in between the elements
e = 12; % number of elements (an even number are required for the even degree cases)
a = 0.0; % left boundary
b = e; % right boundary
% compute the quadrature rule
[nodes, weights] = quadrule(p,r,e,a,b);
% plot the quadrature nodes and weights
stem(nodes,weights,'o')
xlabel('x')
ylabel('w')
title('quadrature nodes versus weights')
% test rule
knots = [a*ones(1,p) linspace(a,b,e+1) b*ones(1,p)];
dim = length(knots)-p-1;
target = zeros(dim,1);
for k=1:dim
target(k) = (knots(k+p+1)-knots(k)) / (p+1);
end
B = spcol(knots,p+1,nodes);
err = norm(target' - weights' * B)/dim;
% output result
sprintf('rule (%d,%d) passed: %d',p,r,err<10e-15)