-
Notifications
You must be signed in to change notification settings - Fork 0
/
TestProject.m
134 lines (115 loc) · 2.66 KB
/
TestProject.m
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
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
% Script to compare orthogonal polynomial projections.
%
% Daniel R. Reynolds
% SMU Mathematics
% Math 4315
% remove all existing variables
clear
% set interval and weight functions
a = -1;
b = 1;
wL = @(x) ones(size(x));
wC = @(x) 1./sqrt(1-x.^2);
% set polynomial degrees for tests
nvals = [5, 10, 15, 20];
% set function that we will interpolate
f = @(x) 1./(1+x.^2);
% set evaluation points for plots
x = linspace(a, b, 2001);
% initialize plots
figure(1) % Chebyshev interpolant
plot(x, f(x), 'DisplayName', 'f(x)')
hold on
figure(2)
hold on
figure(3) % Legendre projection
plot(x, f(x), 'DisplayName', 'f(x)')
hold on
figure(4)
hold on
figure(5) % Chebyshev projection
plot(x, f(x), 'DisplayName', 'f(x)')
hold on
figure(6)
hold on
% loop over polynomial degrees
for n = nvals
fprintf('Testing with n = %i\n', n);
% Chebyshev interpolant
figure(1)
t = cos((2*(0:n)+1)/(2*n+2)*pi);
p = Lagrange(t, f(t), x);
e = abs(f(x)-p);
plot(x, p, 'DisplayName', sprintf('p_{%i}(x), error = %.2e',n,norm(e,inf)))
figure(2)
semilogy(x, e, 'DisplayName', sprintf('|f-p_{%i}|',n))
% Legendre projection
figure(3)
p = zeros(size(x));
for k=0:n
pk = @(x) Legendre(x,k);
c = L2InnerProduct(pk,f,wL,a,b) * (k+1/2);
p = p + c*pk(x);
end
e = abs(f(x)-p);
plot(x, p, 'DisplayName', sprintf('p_{%i}(x), error = %.2e',n,norm(e,inf)))
figure(4)
semilogy(x, e, 'DisplayName', sprintf('|f-p_{%i}|',n))
% Chebyshev projection
figure(5)
p = zeros(size(x));
for k=0:n
pk = @(x) Chebyshev(x,k);
if (k==0)
c = L2InnerProduct(pk,f,wC,a,b) / pi;
else
c = L2InnerProduct(pk,f,wC,a,b) * 2 / pi;
end
p = p + c*pk(x);
end
e = abs(f(x)-p);
plot(x, p, 'DisplayName', sprintf('p_{%i}(x), error = %.2e',n,norm(e,inf)))
figure(6)
semilogy(x, e, 'DisplayName', sprintf('|f-p_{%i}|',n))
end
% finalize plots
figure(1)
hold off
xlabel('x')
ylabel('f(x), p(x)')
legend('Location','Northwest')
title('Chebyshev interpolant')
figure(2)
hold off
set(gca,'YScale','log')
xlabel('x')
ylabel('|f(x)-p(x)|')
legend('Location','Northwest')
title('Chebyshev interpolant error')
figure(3)
hold off
xlabel('x')
ylabel('f(x), p(x)')
legend('Location','Northwest')
title('Legendre projection')
figure(4)
hold off
set(gca,'YScale','log')
xlabel('x')
ylabel('|f(x)-p(x)|')
legend('Location','Northwest')
title('Legendre projection error')
figure(5)
hold off
xlabel('x')
ylabel('f(x), p(x)')
legend('Location','Northwest')
title('Chebyshev projection')
figure(6)
hold off
set(gca,'YScale','log')
xlabel('x')
ylabel('|f(x)-p(x)|')
legend('Location','Northwest')
title('Chebyshev projection error')
% end of script