graphing-functionsMATLABpolynomials
First of all i am new to Matlab, and not sure whether to use polyfit or something else for my problem.
I plotted a graph with the following matlab code:
t=transpose(linspace(-1,1,50))
y=1./(1+25*t.^2)
n=10
A=fliplr(vander(t))
A=A(:,1:n)
x=A\y
u=linspace(-1,1,1000)
g=x(n)*ones(1,1000)
for $i=(n-1):-1:1
g=g.*u+x(i) end
plot(u,g,'-',t,y,'o')
Now my question is: plot the figure with polynomial fit of polynomial degree 2.
How should i do this?
Thanks
With the least squares method you try to 'solve' an system of linear equations $Ax = b$ for $x$, but if $A$ is not square (thats why you cannot solve exactly) the least square approach is $A^tAx = A^tb$. If you solve this straight forward this can generate huge numericl errors (everything is finde though if you solve it exactly, but Matlab can't do that.). That's why one normally uses the qr decomposition $A=QR$ where $Q$ is orthogonal and $R$ is a upper triangular matrix, so you can solve the least square problem by $Ax = QRx = b$ Since $Q$ is orthonormal you can solve this by $x = R^{-1} Q^t b$. I just applied it to your code and it seems to work very well:
Another thing to consider is that polynomials of 'high' degrees tend to oscilate very much even if the points are almost linear. A rule of thumb is that you should expect heavy oscilations f everythiing above degree 7 or 9.
%%
x=[1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5] % experimental x-values
y=[4.3 6.2 10.1 13.5 19.8 22.6 24.7 29.2] % experimental y-values
degree=8; % specify the degree
A = zeros(length(x),degree);
for exponent=0:degree;
for data=1:length(x);
A(data,exponent+1)=x(data).^exponent; % create matrix A
end;
end;
[q,r] = qr(A);
a=r \ q'*y'; % a are the coëfficients of the polynom
a=flipud(a);
fitpolynoom=polyval(a,x);
error=sum((y-fitpolynoom).^2); % calculates the fit-error using the least-squares method
xx = min(x):0.01:max(x);
fitpolynoom=polyval(a,xx);
figure;
plot(x,y,'black*',xx,fitpolynoom,'g-');
error % displays the value of the fit-error in the Matlab command window
First of all, you need to define x0 before you use it. Then you are just plotting a single point each time which is visible but only one pixel. I recommend using
plot(h,err,'o-');
loglog(h,err,'o-');
This shows you where the points are. As @mathreadler pointed out, if you want to see all points of each iteration, you have to plug the plot commands into the loop and use hold, or you save the results for each $i$ in a vector:
f=@(x)cos(x)
x0=pi/4;
f_exact=-sin(x0)
for i=1:21;
h(i)=10^(1-i);
f_approx=(f(x0+h(i))-f(x0-h(i)))./(2*h(i));
err(i)=abs(f_approx-f_exact);
d_error(i)=abs(((h(i)^2)*sin(x0))./6);
end
plot(h,err,'o-');
figure; %opens a new plot window
loglog(h,err,'o-');
The loglog plot loooks decent enough, The problem you are running into (I think) is a loss of accuracy due to cancellation (i.e. when you substract two numbers that are very close).
Best Answer
So you want to fit y as a function of t, right? Use
p = polyfit(t,y,2);
fit = polyval(p,t);
plot(u,g,'-',t,y,'o',t,fit)
The first line is the built-in polynomial fit function. The number 2 is the degree which you specify and it returns the coefficients of the polynomial in p. Note degree 2 means three coefficients. The second line then evaluates the polynomial using the coefficients in p. And then the third, plots them all together.