1> % retta dei minimi quadrati, retta di regressione


1> x = [1 3 4 6]
x =

   1   3   4   6



2> % media delle x


2> mean(x)
ans =  3.5000


3> y = [2 1 3 6]
y =

   2   1   3   6



4> mean(y)
ans =  3


5> sum(y)
ans =  12


6> prod(y)
ans =  36


7> % il punto (3.5,3)  giacce sulla retta di regressione


7> plot(x,y,'ro')


8> plot(x,y,'rx')


9> plot(x,y,'r+')


10> plot(x,y,'ro')


11> p = polyfit(x,y,3);


12> % polinomio di interpolazione di grado al massimo 3


12> hold on


13> xx = linspace(1,6);


14> yp = polyval(p, xx);


15> plot(xx, yp)


16> r = polyfit(x,y,1)
r =

   0.846154   0.038462



17> yr = polyval(r, xx);


18> plot(xx, yr, 'g')


19> % specchio retta di equazione y = x


19> % ci fa scabiare le x con le y


19> % domanda: e' la retta specolare della retta di regressione


19> % la retta di regressione per i dati (y,x)? no!


19> ryx = polyfit(y,x,1)
ryx =

   0.78571   1.14286



20> yrxy = polyval(ryx, xx);


21> plot(xx, yrxy)


22> % plotare la retta bisettrice y = x


22> plot(xx, xx, 'k')


23> hold off


24> % Interpolazione trigonometrica


24> % Quarteroni-Saleri, Esempio 3.6, pagg. 90-91


24> % prendiamo invece un esempio facile, facile


24> xx = 2*pi/100*[0 : 99];


25> f = @(x) 2*cos(3*x) - sin(5*x)
f =

@(x) 2 * cos (3 * x) - sin (5 * x)


26> plot(xx, f(xx))


27> nodi = 2*pi/8*[0 : 7];


28> fnodi = f(nodi);


29> % plotare il polinomio trigonometrico interpolante


29> help interpft
'interpft' is a function from the file C:\Octave\Octave3.6.4_gcc4.6.2\share\octave\3.6.4\m\general\interpft.m

 -- Function File:  interpft (X, N)

 -- Function File:  interpft (X, N, DIM)

     Fourier interpolation.  If X is a vector, then X is resampled with

     N points.  The data in X is assumed to be equispaced.  If X is an

     array, then operate along each column of the array separately.  If

     DIM is specified, then interpolate along the dimension DIM.



     `interpft' assumes that the interpolated function is periodic, and

     so assumptions are made about the endpoints of the interpolation.



     See also: interp1






Additional help for built-in functions and operators is
available in the online version of the manual.  Use the command
'doc <topic>' to search the manual index.

Help and information about Octave is also available on the WWW
at http://www.octave.org and via the help@octave.org
mailing list.


30> z = interpft(fnodi, 100);


31> hold on


32> plot(xx, z, 'r')


33> % trovare c_k della lezione


33> y = fft(nodi); % e' sbagliato, deve essere y = fft(fnodi);  per i valori corretti si veda riga 54 in poi


34> length(y)
ans =  8


35> c = fftshift(y);


36> fftshift([1 2 3 4 5 6 7 8])
ans =

   5   6   7   8   1   2   3   4



37> % a_0 = 2*c_0


37> c(5)
ans =  21.991


38> c = fftshift(y)/8;


39> c(5)
ans =  2.7489


40> % a_4


40> 2*c(1)
ans = -0.78540


41> % a_3, a_2, a_1


41> c(2)+c(8), c(3) + c(7), c(4) + c(6)
ans = -0.78540
ans = -0.78540
ans = -0.78540


42> % b_3, b_2, b_1


42> i*(c(8)-c(2))
ans = -0.32532


43> i*(c(7)-c(3))
ans = -0.78540


44> i*(c(6)-c(4))
ans = -1.8961


45> % Effetto dell'aliasing


45> % Quarteroni-Saleri, Esempio 3.7, pagg. 91-92


45> f=inline('sin(x)+sin(5*x)')
f = f(x) = sin(x)+sin(5*x)


46> g=inline('sin(x)-sin(3*x)')
g = f(x) = sin(x)-sin(3*x)


47> x=0 : 0.01 : 2*pi;


48> xnodi=2*pi/8*[0:7];


49> hold off


50> plot(x,f(x),x,g(x),'g',xnodi,f(xnodi),'o')


51> hold on


52> z=interpft(f(xnodi), 200);


53> plot(2*pi/200*[0:199],z,'r')


54> y = fft(fnodi);


55> c = fftshift(y);


56> 2*c(1)
ans = 4.2188e-015


57> % a_3, a_2, a_1


57> c(2)+c(8), c(3) + c(7), c(4) + c(6)
ans =  16.000
ans = -2.2204e-015
ans = -4.8572e-017


58> c = fftshift(y)/8;


59> c(2)+c(8), c(3) + c(7), c(4) + c(6)
ans =  2.0000
ans = -2.7756e-016
ans = -6.0715e-018


60> % b_3, b_2, b_1


60> i*(c(8)-c(2))
ans =  1.0000


61> i*(c(7)-c(3))
ans = -4.7184e-016


62> i*(c(6)-c(4))
ans = 4.0593e-016


63> quit