% 10.34 - Fall 2006 % HW Set#2, Problem #2a % Rob Ashcraft - Sept. 11, 2006 % This function can do two things: (1) SOlve for the polynomial % coefficients of an Nth-order polynomial given N+1 data points, (2) Find % the best fit coefficients for an Nth-order polynomial given more than N+1 % data points. % f(x) = a0 + a1*x + a2*x^2 + ... + aN*x^N %If you want to always use the largest polynomial possible for the data, %N_order does no have to be passed to the function function a_vec = calc_poly_coeff(x,f,N_order) % Input: % x = data for the independent variable % f = data for dependent variable % N_order = order of the polynomial you want to fit % % Output % a_vec = polynomial coefficients = [a0; a1; ... aN] % % First, check to see if the number of arguments passed to the function is % equal to 2 (i.e. the user only specified the x and f data). If this is % true, the user wants to largest polynomial possible and the coefficients % can be found by solving a linear system: A*x = f. "nargin" in a built-in % matlab command that tell you have many inputs were passed to a function if(nargin == 2) % Generate the system matrix, with rows of the form: [1 x x^2 ... x^N] for i=1:length(x) for j=1:length(x) X_mat(j,i) = x(j)^(i-1); end end % determine and display the condition number of the system matrix disp(['The condition number for the polynomial is: ',num2str(cond(X_mat),'%5.3e')]); disp(' '); % use the backslash command solve for the parameters a_vec = X_mat\f; % If the user supplied the desired order, there are three possibilities: % (1) you have too many parameters for the amount of data, (2) the number % of parameters and data is equal and you can solve a linear system, or (3) % the data set is larger than the number of parameters and a linear % regression must be performed to find the parameter. elseif(nargin == 3) if(N_order+1 > length(f)) disp(['You do not have enough data to solve for the coefficients of the order ',num2str(N_order),' polynomial.']); return; % kicks out of the function... probably would cause errors. elseif(N_order+1 == length(f)) for i=1:length(x) for j=1:length(x) X_mat(j,i) = x(j)^(i-1); end end a_vec = X_mat\f; else for i=1:(N_order+1) for j=1:length(x) X_mat(j,i) = x(j)^(i-1); end end % in this case, a linear regression is performed to get the parameters a_vec = X_mat'*X_mat\X_mat'*f end end