f = @(x) x.^7;
N = 10;
xi = gauss_legendre_points(N+1);
%% Part (a)
D = DerivMatrix(xi,N);
f_dash_xi = D*f(xi);
figure(1)
fplot(@(x)7.*x.^6,[-1,1],’c’)
hold on
plot(xi,f_dash_xi,’xm’)
hold off
xlabel(‘$\xi$’,’Interpreter’,’latex’)
ylabel(‘$df/d\xi$’,’Interpreter’,’latex’)
legend(‘Exact’,’Spectral’)
%% Part (c)
x = (1-xi)+5*(1+xi);
f_dash_x = D*f(x)/4;
figure(2)
fplot(@(x)7.*x.^6,[2,10],’c’)
hold on
plot(x,f_dash_x,’xm’)
hold off
xlabel(‘$x$’,’Interpreter’,’latex’)
ylabel(‘$dg/dx$’,’Interpreter’,’latex’)
legend(‘Exact’,’Spectral’)
function D=DerivMatrix(x,n)
D=zeros(n+1,n+1);
for i=1:n+1
num(i)=1.0;
for k=1:n+1
if k ~= i
num(i)=num(i)*(x(i)-x(k));
end
end
end
for j=1:n+1
den(j)=1.0;
for k=1:n+1
if k ~= j
den(j)=den(j)*(x(j)-x(k));
end
end
end
for i=1:n+1
for j=1:n+1
if i ~= j
D(i,j)=num(i)/(den(j)*(x(i)-x(j)));
end
end
end
for i=1:n+1
for k=1:n+1
if i~=k
D(i,i)=D(i,i)+1./(x(i)-x(k));
end
end
end
end
function x = gauss_legendre_points(k)
syms t
x = double(vpasolve(legendreP(k,t) == 0));
end