NUMERICAL METHODS IN STRUCTURAL
ENGINEERING
(CE60125):COMPUTER ASSIGNMENT 1
Prepared by: Kanishka Bhattacharya
16CE91P01
CIVIL ENGINEERING DEPARTMENT
CODE: POLYNOMIAL APPROXIMATION USING LAGRANGE
INTERPOLATION
MAIN PROGRAMME:
%%------------------------Polynomial
Approximation-----------------------%%
%%---------------------------------Using--------------------------------%%
%%-------------------------Lagrange
Interpolation------------------------%%
%This code performs a polynomial approx. by Lagrange Interpolation of
the
%following functions
% 1-------f(x)=1/(1+x^2)
for -5<x<5
% 2-------f(x)=2
for 0<=x<1
%
f(x)=1
for 1<=x<2
%
f(x)=2
for 2<=x<3
%
f(x)=1
for 3<=x<4
%
f(x)=2
for 4<=x<=5
%The user will be asked for the following
%
a) Which function is to be approximated (1 or 2)
%
b) Number of data points to be used in the approximation
%The ouput will be a plot of the exact function, the approximated
%polynomial and its derivative
c=input('Which function is to be approximated(1 or 2)?\n');
n=input('How many data points are to be used in approximation?\n');
[y,x,xx,yy]=func(c,n);
%Calling subroutine 'func' to return the
%vector of data points (x), vector of
%corresponding function values (y),vector
of
%plotting points (xx) and vector of exact
%functional values (yy)
p=interp_lag(x,y,xx);
%Calling the subroutine interp_lag to
%estimate the value of the approx. poly.
at
%each plotting point (xx)
%
%----------------------------------------------------------------------%%
%%-------------Estimating the derivative of the
approximation------------%%
%
%----------------------------------------------------------------------%%
%Calculating the derivative at each point of xx
%Calculating derivatives at first n-1 points by forward difference
%and the last point by backward difference
k=length(xx);
dp=zeros(size(xx));
%Vector of the 1st deriv. at xx
for i=1:k-1
dp(i)=(p(i+1)-p(i))/(xx(i+1)-xx(i));
%Forward difference
end
dp(k)=(p(k)-p(k-1))/(xx(k)-xx(k-1));
%Backward difference
%
%----------------------------------------------------------------------%%
%%-------------------------------Plots----------------------------------%%
%
%-----------------------------------------------------------------------
%%
%Plot of exact functional values verses approx. functional values
figure
plot(xx,yy,'b')
%Plotting the exact functional values
hold on
plot(xx,p,'r-.')
%Plotting the approximated functional
values
scatter(x,y,'o')
%Plotting the data points
legend('Exact func.','Approx. func.','Data points')
xlabel('x')
ylabel('y')
title('Graph of exact function and approximating
polynomial','FontSize',14)
%Plot of error of the approx. functional values
e=zeros(size(xx));
for i=1:length(xx)
e(i)=abs((p(i)-yy(i)));
end
figure
plot(xx,e,'b')
legend('|fapprox. - fexact|')
xlabel('x')
ylabel('Error')
title('Graph of absolute error','FontSize',14)
%Plot of approx. functional values and its first derivative
figure
plot(xx,p,'r-.')
hold on
plot(xx,dp,'r--','Linewidth',2)
legend('Approx. func.','Deriv. of approx. func.')
xlabel('x')
ylabel('y/dydx')
title('Graph of approximating polynomial and derivative','FontSize',14)
SUBROUTINE:interp_lag
function [yy]=interp_lag(x,y,xx)
% Function to generate lagrangian poly., interpolating the func. values
(y)
% at points (x) to obtain the approximated values (yy) at points (xx)
% Input: x - vector containing the data points
%
y - vector containing the func. values at the data points
%
xx - vector containing the plotting points
% Output: yy - vector containing the approx. func. values at xx
t=length(xx);
yy=zeros(size(xx));
n=length(x);
for i=1:t
s=0;
for j=1:n
p=1;
for k=[1:n]
if k==j
continue
end
p=p*(xx(i)-x(k))/(x(j)-x(k));
end
s=s+p*y(j);
end
yy(i)=s;
end
end
SUBROUTINE:func
function [y,x,xx,yy]=func(c,n)
% Function to generate the data points (x), func. values (y) at data
% points, vector containing the plotting points (xx) and vector
containing
% the approx. func. values at xx
% Input: c - choice of function
%
n - number of data points
% Output: x - vector containing the data points
%
y - vector containing the func. values at the data points
%
xx - vector containing the plotting points
%
yy - vector containing the approx. func. values at xx
if c==1
x=-5:(10/(n-1)):5;
y=zeros(size(x));
for i=1:length(x)
y(i)=1/(1+x(i)^2);
end
d=(5+5)/((100*length(x))-1);
xx=-5:d:5;
yy=zeros(size(xx));
for i=1:length(xx)
yy(i)=1/(1+xx(i)^2);
end
else
x=0:(5/(n-1)):5;
y=zeros(size(x));
for i=1:n
if x(i)>=0 && x(i)<1
y(i)=2;
end
if x(i)>=1 && x(i)<2
y(i)=1;
end
if x(i)>=2 && x(i)<3
y(i)=2;
%for first function
%vector of data points
%vector of function values
%for n data points, 100*n
%plotting points considered
%vector of approx. func. values
%for second function
%vector of data points
%vector of function values
end
if x(i)>=3 && x(i)<4
y(i)=1;
end
if x(i)>=4 && x(i)<=5
y(i)=2;
end
end
d=(5-0)/((100*length(x))-1);
xx=0:d:5;
yy=zeros(size(xx));
for i=1:length(xx)
if xx(i)>=0 && xx(i)<1
yy(i)=2;
end
if xx(i)>=1 && xx(i)<2
yy(i)=1;
end
if xx(i)>=2 && xx(i)<3
yy(i)=2;
end
if xx(i)>=3 && xx(i)<4
yy(i)=1;
end
if xx(i)>=4 && xx(i)<=5
yy(i)=2;
end
end
end
end
%for n data points, 100*n
%plotting points considered
%vector of approx. func. values
for
No. of data points =3
No. of data points =5
No. of data points =7
No. of data points =10
