creating loop that stops at correct formula for the given input.
Show older comments
I am writting a program the finds the polynomial using the least-square method and Gausian elimination and from there displays the correct polynomial equation. The code for the least-square and the Gaussian elimination works fine, where I am having trouble is at the end getting the display of the polynomial equation depending on what value of "n" you show correcty and nothing lower then or higher then the corresponding the chosen value of "n". What loop would for best to do this? My fist thought is and if else loop but how would I set this up to run correctly?
clear, close, clc
%input
x=[0 1 2 3 4 5];
y=[2.1 7.7 13.6 27.2 40.9 61.1];
n=8; %order of the polynomial you want to find
%%%LEAST SQUARE METHOD-FINDS POLYNOMIAL FOR GIVEN DATA SET%%%%%
%INPUT SECTION
x=x; %x-cordinates from data-input; independent vairiables
y=y; %y-cordinates from data-output; dependent vairiables
%CALCULATIONS SECTION
k=length(x); %NUMBER OF AVAILABLE DATA POINTS
m=n+1; %SIZE OF THE COEFFICENT MATRIX
A=zeros(m,m); %COEFFICENT MATRIX
for j=1:m
for i=1:m
A(j,i)=sum(x.^(i+j-2));
end
end
B=zeros(m,1); %FORCING FUNCTION VECTOR
for i=1:m;
B(i)=sum(y.*x.^(i-1));
end
a1=A\B %COEFFICIENTS FOR THE POLYNOMINAL--> y=a0+a1*x+a2*x^2....an*x^n CAN BE REPLACED BY GAUSSIAN ELIMINATION
%%%%%=========GAUSSIAN ELIMINATION TO FIND "a"========%%%%%%
%%%INPUT SECTION
%CALCULATION SECTION
AB=[A B]; %Augumentent matrix
R=size(AB,1); %# OF ROWS IN AB
C=size(AB,2); %# OF COLUMNS IN AB
%%%%FOWARD ELIMINATION SECTION
for J=1:R-1
[M,I]=max(abs(AB(J:R,J))); %M=MAXIMUM VALUE, I=LOCATION OF THE MAXIMUM VALUE IN THE 1ST ROW
temp=AB(J,:);
AB(J,:)=AB(I+(J-1),:);
AB(I+(J-1),:)=temp;
for i=(J+1):R;
if AB(i,J)~=0;
AB(i,:)=AB(i,:)-(AB(i,J)/AB(J,J))*AB(J,:);
end
end
end
%%%%BACKWARDS SUBSTITUTION
a(R)=AB(R,C)/AB(R,R);
for i=R-1:-1:1
a(i)=(AB(i,C)-AB(i,i+1:R)*a(i+1:R)')/AB(i,i);
end
disp(a)
%========END OF GAUSSIAN ELIMINATION=======%%%%%%%%
%STANDARD DEVIATION
Y_bar=mean(y); %ADVERAGE OF y
St=sum((y-Y_bar).^2);
SD=sqrt(St/(k-1)); %STANDARD DEVIATION
%STANDARD ERROR
for i=1:m;
T(:,i)=a(i)*x.^(i-1); %T=INDIVIDUAL POLYNOMIAL TERMS
end
for i=1:k
y_hat(i)=sum(T(i,:));
end
Sr=sum((y-y_hat).^2);
Se=sqrt(Sr/(k-(n+1))); %STANDARD ERROR-Se
%COEFFICIENT OF DETERMINATION
Cd=(St-Sr)/St %COEFFICIENT OF DETERMINATION (r^2)
fprintf('For n=%d. Coefficient of Determination=%0.5f\n',n,Cd)
%EQUATION FOR THE POLYNOMIAL
syms x y
for n=1
a0=a(:,1);a1=a(:,2);a2=a(:,3);
sympref("FloatingPointOutput",true);
y=a0+a1*x+a2*x.^2
end
for n=2
a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);
sympref("FloatingPointOutput",true);
y=a0+a1*x+a2*x.^2+a3*x.^3
end
for n=3
a0=a(:,1);a1=a(:,2);a2=a(:,2);a3=a(:,4);a4=a(:,5);
sympref("FloatingPointOutput",true);
y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4
end
for n=4
a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);
sympref("FloatingPointOutput",true);
y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5
end
for n=5
a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7)
sympref("FloatingPointOutput",true);
y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6
end
for n=6
a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);
sympref("FloatingPointOutput",true);
y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7
end
for n=7
a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);a8=a(:,9);
sympref("FloatingPointOutput",true);
y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7+a8*x.^8
end
for n=8
a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);a8=a(:,9);a9=a(:,10);
sympref("FloatingPointOutput",true);
y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7+a8*x.^8+a9*x.^9
end
for n=9
a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);a8=a(:,9);a9=a(:,10);a10=a(:,11);
sympref("FloatingPointOutput",true);
y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7+a8*x.^8+a9*x.^9+a10*x.^10
end
Accepted Answer
John D'Errico
on 19 Jul 2024
Edited: John D'Errico
on 19 Jul 2024
You have a fundamental problem here, that you fail to understand. You CANNOT solve for more than 6 parameters, given only 6 data points. That means you would be limited to a degree 5 polynomial. 6 or above will yield meaningless garbage.
When you try to use Gaussian elimination for those polynomials of degree higher than 5, the result wil l be a singular system of equations.And even for the degree 5 polynomial, the use of least squares is meaningless, since that will yield an exact interpolation.
Worse, even a degree 5 polynomial may be somewhat suspect, but you will probably survive there.
As far as a loop to determine the "correct" order" Again, you are probably working under a misconception. There is no perfect order, since no polynomial, even an interpolating one, will actually give you zero residuals. So you will need to employ a tolerance, deciding when to stop when the errors are sufficiently small. Even there though, I would not be surprised to see this always require a degree 5 polynomial.
x=[0 1 2 3 4 5];
y=[2.1 7.7 13.6 27.2 40.9 61.1];
Nmax = 5
resnorm = zeros(1,Nmax);
P = cell(Nmax,1);
for n = 1:Nmax
P{n} = polyfit(x,y,n);
resnorm(n) = norm(y - polyval(P{n},x));
end
resnorm
tol = 0.01; % based on 2 significant digits in your data past the decimal point
plot(1:Nmax,resnorm,'-o')
You should see the error in the fit goes to essentially zero only when n==5, where the polynomial fit will be "exact". But even there, it will not be truly exact, due to floating point issues in double precision arithmetic.
resnorm(end)
Nfit = find(resnorm <= tol,1,'first')
xint = linspace(x(1),x(end));
plot(x,y,'ro',xint,polyval(P{Nfit},xint),'-b')
And even that fit is, in my eyes, complete crap, the result of overfitting. The curve passed exactly through the data points, but it also lacks the essential character of that curve, as we can see.
Nfitminus1 = Nfit - 1;
plot(x,y,'ro',xint,polyval(P{Nfitminus1},xint),'-b')
I hope you see and understand the difference.
Anyway, you can do the same fits using whatever tool you want, but you can still use a loop as I show.
More Answers (0)
Categories
Find more on Creating and Concatenating Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
