I'm trying to modify the code to plot the error E for n=1,2,3... I want to use for loop but I'm not sure how to use it. Please help

clear
n = 3; % the order of the polynomial
a = 2.0; % left end of the interval
b = 3.0; % right end of the interval
h = (b - a)/n; % interpolation grid size
t = a:h:b; % interpolation points
f = 1./t; % f(x) = 1./x, This is the function evaluated at interpolation points

%%%% pn(x) = \sum f(t_i)l_i(x)
hh = 0.01; % grid to plot the function both f and p
x = a:hh:b;
fexact = 1./x; %exact function f at x
l = zeros(n+1, length(x)); %%%% l(1,:): l_0(x), ..., l(n+1): l_n(x)
nn = ones(n+1, length(x));
d = ones(n + 1, length(x));

for i = 1:n+1
    for j = 1:length(x)
        nn(i,j) = 1;     
        d(i,j) = 1;     
        for k = 1:n+1
            if i ~= k
                nn(i,j) = nn(i,j) * (x(j) - t(k));
                d(i,j) = d(i,j) * (t(i) - t(k));  
            end      
        end     
        l(i,j) = nn(i,j)/d(i,j);  
    end
end

fapp = zeros(length(x),1);
for j = 1:length(x)
   for i=1:n+1 
       fapp(j) = fapp(j) + f(i)*l(i,j);  
   end
end 
En = 0;
Ed = 0;
for i = 1:length(x) 
     Ed = Ed + fexact(i)^2; 
     En = En + (fexact(i) - fapp(i))^2;
end
Ed = sqrt(Ed);
En = sqrt(En);
E = En/Ed;
display(E)

plot(x,fexact,'b*-')
hold on
plot(x,fapp,'ro-' )

Answer 1

[Hrishikesh Borate]

Hi,

The following code demonstrates the plot of E for n = 1 to 10.

clear
e = [];
figure;
subplot(2,1,1);
ax1 = gca;
subplot(2,1,2);
ax2 = gca;
for n=1:10
    % n = 3; % the order of the polynomial
    a = 2.0; % left end of the interval
    b = 3.0; % right end of the interval
    h = (b - a)/n; % interpolation grid size
    t = a:h:b; % interpolation points
    f = 1./t; % f(x) = 1./x, This is the function evaluated at interpolation points
    
    %%%% pn(x) = \sum f(t_i)l_i(x)
    hh = 0.01; % grid to plot the function both f and p
    x = a:hh:b;
    fexact = 1./x; %exact function f at x
    l = zeros(n+1, length(x)); %%%% l(1,:): l_0(x), ..., l(n+1): l_n(x)
    nn = ones(n+1, length(x));
    d = ones(n + 1, length(x));
    
    for i = 1:n+1
        for j = 1:length(x)
            nn(i,j) = 1;
            d(i,j) = 1;
            for k = 1:n+1
                if i ~= k
                    nn(i,j) = nn(i,j) * (x(j) - t(k));
                    d(i,j) = d(i,j) * (t(i) - t(k));
                end
            end
            l(i,j) = nn(i,j)/d(i,j);
        end
    end
    
    fapp = zeros(length(x),1);
    for j = 1:length(x)
        for i=1:n+1
            fapp(j) = fapp(j) + f(i)*l(i,j);
        end
    end
    En = 0;
    Ed = 0;
    for i = 1:length(x)
        Ed = Ed + fexact(i)^2;
        En = En + (fexact(i) - fapp(i))^2;
    end
    Ed = sqrt(Ed);
    En = sqrt(En);
    E = En/Ed;
    e = [e;E];
    plot(ax1,e);
    
    plot(ax2,x,fexact,'b*-')
    hold on
    plot(ax2,x,fapp,'ro-' )
    hold off
end

[ebtisam almehmadi]

Thank you so much, it works. Could you please explain the steps you do? Iknow little about matlab.