Download numerical analysis for electrical engineering and more Exercises Numerical Methods in Engineering in PDF only on Docsity!
Experiment : Least Squares Regression and Interpolation
Objective : In this lab, students will get acquainted with least squares linear regression and interpolation
techniques.
Introduction : First, we will come up with a Matlab function to compute slope and intercept of the best fit
straight line along with correlation coefficient and standard error of the estimate. For given inputs vectors x
and y, the output slope and intercept along with the best fit plot are provided.
For a given set of n data points (xi, yi), the overall error calculated
The function E has a minimum at the values of a 1 and a 0 where the partial derivatives of E with respect
to each variable is equal to zero. Taking the partial derivatives and setting them equal to zero gives:
Experiment - 1 curve fitting
Exercise Write Matlab Program that best fit Linear regression line for the following data
x 10 20 30 40 50 60 70 80
y 25 70 380 550 610 1220 830 1450
%Linear regression of y=a1x+a clc; close all; x=[10 20 30 40 50 50 70 80]; y=[25 70 380 550 610 1220 830 1450]; linreg(x,y); function[a1,a0]=linreg(x,y) n=length(x); sx=sum(x); sy=sum(y); sxx=sum(x.x); sxy=sum(x.y); a1=(n.sxy-sx.sy)./(n.sxx-sx.^2); a0=(sy-a1.sx)./n; xx=linspace(min(x),max(x),100); yy=a1.xx+a0; plot(x,y,'ko',xx,yy,'r-') grid; syms x; y=a1.x+a end
Exercise Implement the above problems using built in MATLAB function polyfit,
clc; close all; x=[10 20 30 40 50 50 70 80]; y=[25 70 380 550 610 1220 830 1450]; a=polyfit(x,y,1)
Exercise The following data is given:
y = (12159*x)/638 - 61220/
outputs
a =19.0580 - 191.
Exercise The set of the following five data points is given:
Develop a MATLAB program that finds the interpolation of y at x=3 in Lagrange interpolation polynomial
%Lagrange interpolation polynoial clc; x=[1 2 4 5 7]; y=[52 5 - 5 - 40 10]; x0=3; LagrangeINT(x,y,x0); function Yint = LagrangeINT(x,y,Xint) n = length(x); for i =1:n L(i)=1; for j=1:n if j~=i L(i)= L(i)(Xint-x(j))/(x(i)-x(j)); end end end Yint = sum(y.L) end
Exercise The following data is given:
Write the polynomial in Lagrange form that passes through the points; then use it to calculate the
interpolated value of y at x = 5..
Output
Yint =6.
Exercise find the interpolation of the newton divided difference at x=0.5 for the given data points
clc; x1 = [- 2 - 1 1 2 4]; y1 = [- 6 0 0 6 60]; x0=0.5; interpolated=Newtint(x1,y1,x0) function yint = Newtint(x,y,xx) % Newtint: Newton interpolating polynomial % yint = Newtint(x,y,xx): Uses an (n - 1)-order Newton % interpolating polynomial based on n data points (x, y) % to determine a value of the dependent variable (yint) % at a given value of the independent variable, xx. % input: % x = independent variable % y = dependent variable % xx = value of independent variable at which % interpolation is calculated % output: % yint = interpolated value of dependent variable % compute the finite divided differences in the form of a % difference table n = length(x); if length(y)~=n, error('x and y must be same length'); end b = zeros(n,n); % assign dependent variables to the first column of b. b(:,1) = y(:); % the (:) ensures that y is a column vector. for j = 2:n for i = 1:n-j+ b(i,j) = (b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i)); end end disp('Newton divided difference table'); disp(b); disp('Newton Coefficents are'); disp(b(1,:)); % use the finite divided differences to interpolate xt = 1; yint = b(1,1); for j = 1:n- 1 xt = xt(xx-x(j)); yint = yint+b(1,j+1)xt; end N=length(x)-1; a = b(1,:); m = a(N+1); %Begin for k = N:-1: m = [m a(k)] - [0 mx(k)]; %n(x)(x - x(k - 1))+a_k - 1 end disp('coffients of x in decreasing order of power x') m xx = (-2:0.02: 4 ); yy = polyval(m,xx); plot(xx,yy,'b-',x,y,'*') end
Exercise use polyfit function to check the results you obtained in the above problem {(−2, −6), (−1, 0), (1, 0), (2, 6), (4, 60)} clc; x1 = [- 2 - 1 1 2 4]; y1 = [- 6 0 0 6 60]; a=polyfit(x1,y1,4) xx=(-2:0.02:4); yy=polyval(a,xx); plot(xx,yy,'b-',x1,y1,'*') grid
outputs
Newton divided difference table
Newton Coefficents are
coffients of x in decreasing order of power x
m =
interpolated =
Outputs
a =
- 0.0000 1.0000 0.0000 - 1.0000 - 0.
Table Lookup
. First, it would search the independent variable vector to find the interval containing the unknown. Then it would perform the linear interpolation using one of the techniques described For ordered data, there are two simple ways to find the interval ➢ sequential search this method involves comparing the desired value with each element of the vector in sequence until the interval is located ➢ binary search. For situations for which there are lots of data, the sequential sort is inefficient because it must search through all the preceding points to find values The approach is akin to the bisection method for root location. Just as in bisection, the index at the midpoint iM is computed as the average of the first or “lower” index iL = 1 and the last or “upper” index iU = n. The unknown xx is then compared with the value of x at the midpoint x(iM) to assess whether it is in the lower half of the array or in the upper half. Depending on where it lies, either the lower or upper index is redefined as being the middle index. The process is repeated until the difference between the upper and the lower index is less than or equal to zero. At this point, the lower index lies at the lower bound of the interval containing xx, the loop terminates, and the linear interpolation is performed
Exercise use both sequential and binary search methods to find the linear spline interpolation of the following data at T= T = [-40 0 20 50 100 150 200 250 300 400 500]; density = [1.52 1.29 1.2 1.09 .946 .935 .746 .675 .616 .525 .457]; clc T = [- 40 0 20 50 100 150 200 250 300 400 500]; density = [1.52 1.29 1.2 1.09 .946 .935 .746 .675 .616 .525 .457]; TableLookBin(T,density,350)% binary search TableLook(T,density,350)% sequential search function yi = TableLookBin(x, y, xx) n = length(x); if xx < x(1) | xx > x(n) error('Interpolation outside range') end % binary search iL = 1; iU = n; while (1) if iU - iL <= 1, break, end iM = fix((iL + iU) / 2); if x(iM) < xx iL = iM; else iU = iM; end end % linear interpolation yi = y(iL) + (y(iL+1)-y(iL))/(x(iL+1)-x(iL))(xx - x(iL)); end function yi = TableLook(x, y, xx) n = length(x); if xx < x(1) | xx > x(n) error('Interpolation outside range') end % sequential search i = 1; while(1) if xx <= x(i + 1), break, end i = i + 1; end % linear interpolation yi = y(i) + (y(i+1)-y(i))/(x(i+1)-x(i))(xx-x(i)); end
Outputs
ans =
ans =
Exercise implement a spline interpolation ( natural, clamped ) for the data given below and find interpolation at x=1. x=[0 1 2 3 4 5] y=[3 1 4 1 2 0] clc xd=[0 1 2 3 4 5]; yd=[3 1 4 1 2 0]; xx=1.5; nd=length(xd); [a,b,c,d]=genspline(xd,yd,nd); evalspline(a,b,c,d,xx,xd,nd) plotspline(a,b,c,d,xd,nd) function[a,b,c,d]=genspline(x,y,n) for i=1:n- 1 h(i)=x(i+1)-x(i); end A(1,1)=1; A(n,n)=1; f=zeros(n,1); for i=2:n- 1 A(i,i)=2(h(i)-h(i-1)); f(i)=6.(y(i+1)-y(i))./h(i)-(y(i)-y(i-1))./h(i-1); end for i=2:n- 2 A(i,i+1)=h(i+1); end for i=3:n- 1 A(i,i-1)=h(i); end s=A\f; for i=1:n- 1 a(i)=(s(i+1)-s(i))/(6h(i)); b(i)=s(i)/2; c(i)=(y(i+1)-y(i))/h(i)-(2h(i)s(i)+h(i)s(i+1))/6; d(i)=y(i); end [a' b' c' d'] end function yy=evalspline(a,b,c,d,xx,x,n) i=1; while(xx>x(i+1)&i<n-1) i=i+1; end yy=a(i).(xx-x(i)).^3+b(i).(xx-x(i)).^2+c(i).(xx-x(i))+d(i); end function plotspline(a,b,c,d,x,n) dx=(x(n)-x(1))/100; for j=1: xx(j)=x(1)+dx(j-1); yy(j)=evalspline(a,b,c,d,xx(j),x,n); end xd=[0 1 2 3 4 5]; yd=[3 1 4 1 2 0]; plot(xx,yy,'-',xd,yd,'*') end
xa = zeros(1,n+1); xa(1) = 3.0 * (a(2) - a(1)) / h(1) - 3.0 * fp0; xa(n+1) = 3.0 * fpn - 3.0 * (a(n+1) - a(n)) / h(n); for i = 1:m xa(i+1) = 3.0(a(i+2)h(i)-a(i+1)(x(i+2)-x(i))+a(i)h(i+1))/(h(i+1)h(i)); end xl = zeros(1,n+1); xu = zeros(1,n+1); xz = zeros(1,n+1); xl(1) = 2.0 * h(1); xu(1) = 0.5; xz(1) = xa(1) / xl(1); for i = 1:m xl(i+1) = 2.0 * (x(i+2) - x(i)) - h(i) * xu(i); xu(i+1) = h(i+1) / xl(i+1); xz(i+1) = (xa(i+1) - h(i) * xz(i)) / xl(i+1); end xl(n+1) = h(n) * (2.0 - xu(n)); xz(n+1) = (xa(n+1) - h(n) * xz(n)) / xl(n+1); c = zeros(1,n+1); b = zeros(1,n+1); d = zeros(1,n+1); c(n+1) = xz(n+1); for i = 1:n j = n - i; c(j+1) = xz(j+1) - xu(j+1) * c(j+2); b(j+1) = (a(j+2)-a(j+1))/h(j+1)-h(j+1)(c(j+2)+2.0*c(j+1))/3.0; d(j+1) = (c(j+2) - c(j+1)) / (3.0 * h(j+1)); end fprintf('The numbers x(0), ..., x(n) are:\n'); for i = 0:n fprintf(' %5.4f', x(i+1)); end fprintf('\n\nThe coefficients of the spline on the subintervals are: \n'); fprintf(' a(i) b(i) c(i) d(i)\n'); for i = 0:m fprintf('%11.8f %11.8f %11.8f %11.8f\n',a(i+1),b(i+1),c(i+1),d(i+1)); end
Exercise use built in libraries to do the problems given above xd=[0 1 2 3 4 5]; yd=[3 1 4 1 2 0]; x=0:0.1:5; yi=interp1(xd,yd,x,'spline'); plot(xd,yd,'*',x,yi,'-'); Exercise :Analyse the out come for different methods of interpolation for the above problems ( nearest, linear,pchip)
Enter f'(x(0)) and f'(x(n)) on separate lines
The numbers x(0), ..., x(n) are:
The coefficients of the spline on the subintervals are:
a(i) b(i) c(i) d(i)
