creating loop that stops at correct formula for the given input.

1 view (last 30 days)
Jon
Jon on 19 Jul 2024
Edited: John D'Errico on 19 Jul 2024
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

Sign in to answer this question.

Accepted Answer

John D'Errico
John D'Errico on 19 Jul 2024
Edited: John D'Errico on 19 Jul 2024
Ran in:
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
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
resnorm = 1x5
11.5327 1.9356 1.8366 1.8268 0.0000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
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)
ans = 1.8375e-13
Nfit = find(resnorm <= tol,1,'first')
Nfit = 5
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.
  1 Comment
Jon
Jon on 19 Jul 2024
I understand that, I am trying to set this up to use on another problem that has an 8 by 8 matrix. I was hoping there was a way to stop it at the correct polynomial based on the number of "n" and the number of data points. Even when I put in the 8 data points it doen't stop at the value of "n" chosen. Any suggestions to allow this to work if the data points and value of "n" create a conflict and if they don't conflict to stop at the correct polynomial? Below is the ammend code with 8 data points to show it still doen't stop where the "n" input value is selected.
clear, close, clc
%input
x=[0 1 2 3 4 5 6 7 8];
y=[2.1 7.7 13.6 27.2 40.9 61.1 65.9 73.6 80.1];
n=input('order of polynomial (1,2,3,4,5,6,7,8,9):'); %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
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.576099e-16.
a1 = 8×1
2.132944824015855
-6.816967090378475
27.072915338895836
-22.177764982376821
8.793838303081925
-1.725555425180783
0.162962951294453
-0.005932539270024
%%%%%=========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)
2.132944824061262 -6.816967090394376 27.072915338897619 -22.177764982376868 8.793838303081923 -1.725555425180783 0.162962951294453 -0.005932539270024
%========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)
Cd =
0.998053904721814
fprintf('For n=%d. Coefficient of Determination=%0.5f\n',n,Cd)
For n=7. Coefficient of Determination=0.99805
%EQUATION FOR THE POLYNOMIAL
syms x y
for i=1:n:9
error=0;
switch n
case 1
a0=a(:,1);a1=a(:,2);a2=a(:,3);
sympref("FloatingPointOutput",true);
y=a0+a1*x+a2*x.^2
case 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
case 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
case 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
case 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
case 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
case 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
case 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
case 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
otherwise
error=1
end
end
Index in position 2 exceeds array bounds. Index must not exceed 8.
if error
disp('error current or polynomial to high producing inaccuries')
else
fprintf('y=',y)
end

Sign in to comment.

More Answers (0)

Categories

Find more on Polynomials in Help Center and File Exchange

Tags

Products


Release

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!