% Class 6 - Interpolation (polynomials) % Coded by NFR on 9.12.2019 % % Problem 17.11 as an example close all clc T = (200:50:450); D = [1.708 1.367 1.139 0.967 0.854 0.759]; % GOAL: Density at 330K % Using a first order polynomial % Two points are needed for 1st order % Pick the two closest to the desired value % Points 3 and 4 p = polyfit(T(3:4),D(3:4),1); %<-- n = order of polynomial % p(1) = slope, p(2) = intercept % Use polyval to find solution Sol1 = polyval(p,330) % Third order fit % 4 points p = polyfit(T(2:5),D(2:5),3); Sol2 = polyval(p,330) % 5th Order fit p = polyfit(T,D,1); Sol3 = polyval(p,330) % Plot all the points plot(T,D,':o',330,Sol1,'r+',330,Sol2,'bo',330,Sol3,'k+') % One more example, to show wrinkle of % inverse interpolation % (this will also show poorly conditioned error) Y_vec = (1920:10:2000); P_vec = [106.46 123.08 132.12 152.27 180.67... 205.05 227.23 249.46 281.42]; plot(Y_vec,P_vec,'o') xlabel('Year') ylabel('Population (Millions)') % Show plot of switching axes, % why you DON'T do this for inverse interpolation figure plot(P_vec,Y_vec,'o') ylabel('Year') xlabel('Population (Millions)') % You can't use polyfit/polyval in reverse % (usually) because the y axis data is not % at even intervals % % Regular interpolation % What is the pop in 1936? % Test an 8th order polynomial p = polyfit(Y_vec(1:4),P_vec(1:4),3); Sol = polyval(p,1936); plot(Y_vec,P_vec,'o',1936,Sol,'r+') xlabel('Year') ylabel('Population') % Inverse interpolation % Conditional logic method to solve % Can use root finding (later) % In what year did the pop hit 190 M? p = polyfit(Y_vec(4:7),P_vec(4:7),3); x_vec = linspace(1950,1980,1000); sol_vec = polyval(p,x_vec); [~,ind] = min(abs(sol_vec-190)); Year = x_vec(ind); fprintf('The year is %.2f\n',Year)