diary('aproximacio.txt') diary on %% 2. Aproximació %% %% 1. Useu els pesos moleculars dels sis òxids de nitrogen donats aquí %% per a calcular els pesos atòmics del nitrogen i l'oxigen %% %% NO : 30.006 N2 O3 : 76.012 N2 O : 44.013 %% N2 O5 : 108.010 NO2 : 46.006 N2 O4 : 92.011 %% %% Solució: >> format long >> A=[1 3 1 5 2 4; 1 2 2 2 1 2]' A = 1 1 3 2 1 2 5 2 2 1 4 2 >> b=[30.006,76.012,44.013,108.010,46.006,92.011]' b = 30.006 76.012 44.013 108.010 46.006 92.011 >> M=A'*A; f=A'*b; >> det(M) ans = 167.000000000000 % ~= 0 >> x=M\f x = 15.9992934131737 14.0069161676647 % Per tant: % Pes atòmic de (O) = 15.999 % Pes atòmic del (N) = 14.007 % D'altra banda, la desviació quadràtica, sigma, és: >> sigma = norm(A*x-b)^2 ans = 4.79041916165687e-07 % Nota important: com que rang(A)= 2, ja que, per exemple, >> det(A(1:2,1:2)) ans = -1 % llavors, A\b dóna directament la solució que minimitza la norma % quadràtica. En efecte: >> x = A\b x = 15.9992934131737 14.0069161676647 % Recordem, si el sistema Ax=f és compatible determinat, la comanda de % MATLAB/Octave x=A\f, retorna la (única) solució. Si és compatible % indeterminat, ens dóna una solució qualsevol i si, com en aquest cas, % és incompatible, torna la solució per mínims quadrats (i.e., el valor % de x que fa mínim norm(Ax-f,2)) %% ----------------------------------------------------------------- %% %% 2. Hom suposa que el cometa Tentax, descobert a l'any 1968, és un %% objecte del Sistema Solar. En un cert sistema de coordenades %% polars (r, phi), centrat en el Sol, s'han mesurat experimental- %% ment les següents posicions del cometa: %% %% r | 2.70 2.00 1.61 1.20 1.02 (UA) %% ----------------------------------------- %% phi | 48 67 83 108 126 (deg) %% %% per les lleis de Kepler i menystinguent les pertorbacions dels %% planetes, el cometa es mou en una cònica que, en les coordenades %% polars usades, té per equació r = p /(1 + e cos(phi)) on p és un %% paràmetre i e és l'excentricitat de la trajectòria. Ajusteu per %% mínims quadrats els valors de p i e a partir de les mesures %% fetes. %% %% Solució: >> format long >> phi=[48 67 83 108 126]; >> r=[2.70 2.00 1.61 1.20 1.02]; % Posem: -e*r*cos(phi) + p = r i definim x = [e, p]'. Aleshores, amb % les dades de la taula, el corresponent sistema d'equacions és de la % forma Ax = r, amb >> A=[-r.*cosd(phi),ones(5,1)] A = -1.806652637168917 1.000000000000000 -0.781462256978547 1.000000000000000 -0.196209642882287 1.000000000000000 0.370820393249937 1.000000000000000 0.599540957338322 1.000000000000000 % D'on les equacions normals són M*x = f, on >> M=A'*A M = 4.41013235800759 -1.81396318644149 -1.81396318644149 5.00000000000000 >> f=A'*r f = -5.70026791096864 8.53000000000000 >> det(M) ans = 18.7601993482730 % ~= 0 >> x=M\f x = -0.694461361131358 1.454054531300334 % Per tant: e = -0.694 % p = 1.45 % Nota: (1) com al problema 1, rang(A)=2. En efecte, >> det(A(1:2,1:2)) ans = -1.02519038019037 % ~= 0 % Llavors la solució es pot obtenir directament, >> A\r ans = -0.694461361131358 1.454054531300333 % (2) També podem fer servir la funció polyfit, que ---segons diu el % guió de les pràctiques---: "de fet és la comanda més versàtil per % resoldre problemes d'aproximació amb MATLAB/Octave". Els detalls del % seu ús els podem trobar al help del MATLAB/Octave: % >> help polyfit % ..... % -- Function File: [P, S, MU] = polyfit (X, Y, N) % Return the coefficients of a polynomial P(X) of degree N that % minimizes the least-squares-error of the fit. % % The polynomial coefficients are returned in a row vector. % ..... >> x=polyfit(-r.*cosd(phi),r,1) x = -0.694461361131357 1.454054531300333 % = e, p % Finalment, dibuixem l'òrbita i les posicions del cometa donades a la % taula >> e=x(1);p=x(2); >> theta=pi/4:0.01:7*pi/4; >> rr=p./(1+e*cos(theta)); >> polar(phi*pi/180,r,'or'); >> hold on;polar(theta,rr,'-b'); >> hold off %% ----------------------------------------------------------------- %% %% 3. El nivell de l'aigua del mar del Nord està determinat principal- %% ment per l'anomenada marea M2, el període de la qual és al voltant %% de les 12 hores i, per tant, suposem que té aproximadament la %% forma: %% %% H(t) = h0 + a1 * sin (2*t*pi/12) + a2 * cos(2*t*pi/12), %% %% on t es mesura en hores. S'han fet les següents mesures: %% %% t | 0 2 4 6 8 10 hores %% --------------------------------------------------- %% H(t) | 1.0 1.6 1.4 0.6 0.2 0.8 metres %% %% (a) Són ortogonals 1, cos(2*pi*t/12), sin (2*pi*t/12)? %% (b) Ajusteu H(t) a aquestes mesures usant mínims quadrats. %% (c) Trobeu la desviació quadràtica mitjana. %% %% Solució: >> format long >> t=[0 2 4 6 8 10]; >> H=[1.0 1.6 1.4 0.6 0.2 0.8]; % (a) Són ortogonals. En efecte: >> format short; A=[ones(1,6); sin(2*pi*t/12); cos(2*pi*t/12)] A = 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 0.00000 0.86603 0.86603 0.00000 -0.86603 -0.86603 1.00000 0.50000 -0.50000 -1.00000 -0.50000 0.50000 >> format long; M=A*A' M = 6.00000000000000e+00 4.44089209850063e-16 -1.11022302462516e-16 4.44089209850063e-16 3.00000000000000e+00 2.92436430973653e-16 -1.11022302462516e-16 2.92436430973653e-16 3.00000000000000e+00 % (b) Ajustem H(t) a les mesures de taula la usant mínims quadrats: >> f=A*H' f = 5.600000000000000 1.732050807568877 0.800000000000001 % Com que M és una matriu diagonal i clarament el seu determinant és % ~= 0, el sistema M x = f és compatible determinat i té per solució: >> x=M\f x = 0.933333333333333 % ~ 5'6/6 0.577350269189626 % ~ sqrt(3)/3 0.266666666666667 % ~ 4/15 >> a0=x(1),a1=x(2),a2=x(3) a0 = 0.933333333333333 a1 = 0.577350269189626 a2 = 0.266666666666667 % (c) Desviació quadràtica, sigma: >> HH=inline('a0 + a1*sin(2*pi*t/12) + a2*cos(2*pi*t/12)','t') HH = f(t) = a0 + a1*sin(2*pi*t/12) + a2*cos(2*pi*t/12) >> sigma=norm(HH(t)-H,2)^2 sigma = 0.120000000000000 % Per últim, dibuixem la funció aproximada (i. e., el polinomi tri- % gonomètric) v.s. els valors de la taula: >> x=0:0.1:10; >> plot(t,H,'or',x,HH(x),'-b') >> hold on >> xlabel('t (h)') >> ylabel('Nivell de l''aigua (m)') >> title('Marea M2') >> hold off %% ----------------------------------------------------------------- %% %% 4. Trobeu la millor aproximació per mínims quadrats als punts %% ( -pi/2 ,1), (0, 0), ( pi/2 , 1 /2), (pi, 1) i que sigui del tipus %% f(theta)= a + r*sin(theta+alpha), amb a, r reals i alpha entre 0 i %% 2*pi. %% %% Solució: % f(theta)= a + r*sin(theta+alpha) % = a + r*sin(alpha)*cos(theta)+r*cos(alpha)*sin(theta) >> format long >> theta=[-pi/2,0,pi/2,pi]; >> b=[1,0,1/2,1]; >> format short,A=[ones(1,4);cos(theta);sin(theta)] A = 1.00000 1.00000 1.00000 1.00000 0.00000 1.00000 0.00000 -1.00000 -1.00000 0.00000 1.00000 0.00000 >> format long; M=A*A' M = 4.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 2.000000000000000 -0.000000000000000 0.000000000000000 -0.000000000000000 2.000000000000000 >> f=A*b' f = 2.500000000000000 -1.000000000000000 -0.500000000000000 % De nou, M és una matriu diagonal i el seu determinant és ~= 0, el % sistema M x = f és compatible determinat i té per solució: >> x=M\f x = 0.625000000000000 -0.500000000000000 -0.250000000000000 >> a0=x(1),a1=x(2),b1=x(3) a0 = 0.625000000000000 a1 = -0.500000000000000 b1 = -0.250000000000000 % Per últim, tenim que: a0=a, a1=r*sin(apha) i b1=r*cos(alpha). D'on a, % r i alpha resulten: >> a=a0,r=sqrt(a1^2+b1^2),alpha=atan2(a1,b1) a = 0.625000000000000 r = 0.559016994374947 alpha = -2.03444393579570 % ...o bé, si volem alpha entre 0 i 2*pi >> alpha=2*pi+alpha alpha = 4.24874137138388 %% ----------------------------------------------------------------- %% %% 5. Per als set planetes més propers al Sol i usant mesures de Tycho %% Brahe, Kepler va obtenir una taula semblant a la següent: %% %% %% r | 60 110 150 230 780 1430 2870 x 1.E6 Km %% ---------------------------------------------------------------- %% T | 90 225 365 690 4330 10750 30650 dies %% %% on r mesura la distancia mitjana del planeta al sol i T el període %% orbital. %% %% Useu aquests valors per a fer dibuixos de gràfiques (r, T ), %% (log(r), T ), etc...fins que trobeu que els resultats es disposin %% aproximadament en una recta. En cas de que sigui així, useu la %% taula per fer el "millor ajust" i determinar la fórmula que %% relacioni r i T. %% %% Solució: >> format long >> r=[60 110 150 230 780 1430 2870]'; >> T=[90 225 365 690 4330 10750 30650]'; % Provem primer amb (r,T): >> plot(r,T,'o-g',r,T,'-b'); xlabel('r');ylabel('T'); % ...no sembla una recta, % a continuació amb (log(r), T), >> plot(log(r),T,'o-r',log(r),T,'-b'); >>xlabel('log r');ylabel('T'); % pitjor! % mirem si fent (log(r), log(T))... >> plot(log(r),log(T),'o-r',log(r),log(T),'-b'); >> xlabel('log r');ylabel('log T'); % Veiem que aquesta última opció és la més acceptable. Així tractarem % d'ajustar per mínims quadrats els valors de la taula per una recta del % tipus % % log(T)= alpha * log(r) + beta % % Podem plantejar les equacions normals i resoldre-les: >> A=[log(r),ones(7,1)],M=A'*A,detM=det(M) A = 4.09434456222210 1.00000000000000 4.70048036579242 1.00000000000000 5.01063529409626 1.00000000000000 5.43807930892320 1.00000000000000 6.65929391968364 1.00000000000000 7.26542972325395 1.00000000000000 7.96206730875367 1.00000000000000 M = 254.06452608538066 41.13033048272523 41.13033048272523 7.00000000000000 detM = 86.7475969794682 % ~= 0 b=log(T),f=A'*b,x=M\f b = 4.49980967033027 5.41610040220442 5.89989735358249 6.53669159759130 8.37332282099653 9.28266103355581 10.33038794136436 f = 314.4455093397307 50.3388708196252 x = 1.50625696837047 % = alpha -1.65913944019430 % = beta % i agafem alpha = 1.506, beta = -1.659 com a valors aproximats. % % Alternativament, fent ús de la funció polyfit del MATLAB/Octave, % trobem: >> x=polyfit(log(r),log(T),1) x = 1.50625696837047 -1.65913944019427 % = alpha, beta. % Més encara, com que rang(A)=2, ja que per exemple, >> detA=det(A([1,4],[1,2])) detA = -1.34373474670110 % ~= 0 % l'observació que fèiem al final del problema 1, també s'aplica % aquí. Així, directament: >> x=A\b x = 1.50625696837047 % = alpha -1.65913944019427 % = beta % Per últim, la desviació quadràtica, sigma, d'aquesta solució és: >> sigma=norm(A*x-b,2)^2 sigma = 2.68744304738384e-04 % De fet, la tercera Llei de Kepler (Llei dels períodes, 1618) estableix % que T^2/r^3 = constant. Això coincideix força bé amb el resultat % obtingut si arrodonim alpha = 3/2, ja que % % log(T)= alpha*log(r)+beta <=> T = exp(beta)*r^alpha % % Finalment, representem els punts de la taula i aquesta última funció, % T=exp(beta)*r^alpha, pels valors de alpha i beta obtinguts. % >>rho=60:0.01:2870; >> P=inline('exp(beta)* r.^alpha','r') P = f(r) = exp(beta)* r.^alpha >> plot(r,T,'or',rho,P(rho),'-b');xlabel('r');ylabel('T'),title(... '3a. Llei de Kepler: T^2/r^3 = constant'); %% ----------------------------------------------------------------- %% %% 6. a) Donades les abscisses x0 = -2 , x1 = -1, x2 = 0, x3 = 1, %% x4 = 2, determineu a, b i c per tal que els polinomis %% %% Psi_0(x) = 1, Psi_1(x) = x + a, Psi_2(x) = x^2 + b x + c, %% %% siguin ortogonals respecte d'elles. %% %% b) Amb els polinomis Psi_i(x) obtinguts, trobeu c0, c1 i c2 %% de manera que, %% %% p(x) = c0 * Psi_0(x) + c1 * Psi_1(x) + c2 * Psi_2(x), %% %% sigui la millor aproximació per mínims quadrats de la taula %% %% %% x | -2 -1 0 1 2 %% --------------------------- %% f(x) | 0 0 3 4 0 %% %% Solució: >> format long % (a) definim el producte escalar discret de dos polinomis p i q % respecte d'un conjunt d'abscisses x=[x(1),...,x(n)], i.e., % n % ----/ % prod = \ p(x(j))*q(x(j)) % / % ----\ % j=1 % % % com una funció inline, ddot >> ddot=inline('sum(polyval(p,x).*polyval(q,x))','p','q','x'); % de manera que si p i q són dos polinomis, aleshores, calcularíem el % producte escalar discret, prod, de p i q respecte dels punts en x % amb la comanda % >> prod=ddot(p,q,x) >> x=[-2 -1 0 1 2] x = -2 -1 0 1 2 % Phi0, Phi1, Phi2: Polinomis de la base canònica >> Phi0=[0 0 1]% Phi0(x)=1 Phi0 = 0 0 1 >> Phi1=[0 1 0]% Phi1(x)=x Phi1 = 0 1 0 >> Phi2=[1 0 0]% Phi2(x)=x^2 Phi2 = 1 0 0 % Ortogonalització de Gram-Schmidt >> Psi0=Phi0 Psi0 = 0 0 1 >> Psi1=Phi1-ddot(Phi1,Psi0,x)/ddot(Psi0,Psi0,x)*Psi0 Psi1 = 0 1 0 >> Psi2=Phi2-ddot(Phi2,Psi0,x)/ddot(Psi0,Psi0,x)*Psi0-... ddot(Phi2,Psi1,x)/ddot(Psi1,Psi1,x)*Psi1 Psi2 = 1 0 -2 % Per tant, tenim: % Psi0 = [0,0,1] --> Psi0(x)=1 % Psi1 = [0,1,0] --> Psi1(x)=x % Psi2 = [1,0,-2] --> Psi2(x)=x^2-2 % (b) Busquem c0, c1, c2 de manera que p(x)=c0*Psi0(x) + % c1*Psi1(x)+c2*Psi2(x) sigui la millor aproximació per mínims % quadrats de la taula >> format rat >> f=[0 0 3 4 0]'; >> A=[polyval(Psi0,x') polyval(Psi1,x') polyval(Psi2,x')] A = 1 -2 2 1 -1 -1 1 0 -2 1 1 -1 1 2 2 >> M=A'*A M = 5 0 0 0 10 0 0 0 14 >> b=A'*f b = 7 4 -10 >> c=inv(M)*b c = 7/5 2/5 -5/7 %% ----------------------------------------------------------------- %% %% 7. La concentració a la sang d'un medicament, mesurat a intervals %% de 30 min ve donada per la taula següent: %% %% temps (min.) | 30 60 90 120 %% ---------------------------------------------- %% Concentració (mg/L) | 60.6 36.8 22.3 13.6 %% %% (a) Assumint que la concentració d'aquest medicament en sang en %% un instant t es pot modelar per un procés de decaïment %% exponencial, %% %% C_t = C_0 * exp(-K*t) %% %% ajusteu els paràmetres C_0 i K per mínims quadrats. %% %% (b) Sabent que la vida mitjana d'una substància es defineix com %% el temps necessari per tal que la seva concentració %% disminueixi a la meitat, calculeu la vida mitjana d'aquesta %% substància. %% %% (c) Calculeu una base de polinomis ortogonals fins a grau 2 %% usant els temps donats a la taula. %% %% Solució: >> format long octave> t=[30 60 90 120] t = 30 60 90 120 >> C_t=[60.6 36.8 22.3 13.6] C_t = 60.6000000000000 36.8000000000000 22.3000000000000 13.6000000000000 % (a) Ajustarem per mínims quadrats les dades a la recta % % log(C_t) = c(1) + c(2)*t, % % amb c(1)=log(C_0), c(2)=-K % >> A=[ones(1,4);t]' A = 1 30 1 60 1 90 1 120 >> M=A'*A M = 4 300 300 27000 >> detM=det(M) detM = 18000.0000000000 >> b=A'*log(C_t)' b = 13.4244492094582 932.0798936937385 >> c=inv(M)*b c = 4.6020089192917091 -0.0166119548923620 % Finalment, tenim C_0 i K >> C_0 = exp(c(1)) C_0 = 99.6843724845930 >> K = -c(2) K = 0.0166119548923620 % Llavors les constants buscades són: % C_0 = 99.6843724845930 mg/L % K = 0.0166119548923620 min^(-1) % A continuació, tot i que l'enunciat no ens ho demana, dibuixem els % punts de la taula i la gràfica de la funció C_t = C_0*exp(-K*t), amb % els valors de K i C_0 trobats >> tt = 30:0.01:120; % xarxa amb 9001 punts >> C_tt = C_0*exp(-K*tt); >> xmin=30; >> xmax=120; >> ymin=min(C_tt); >> ymax=max(C_tt); >> deltaX=5 deltaX = 5 >> deltaY=5 deltaY = 5 >> plot(t,C_t,'or',tt,C_tt,'b-') >> axis([xmin,xmax,ymin-deltaY,ymax+deltaY]) >> hold on >> xlabel('t') >> ylabel('C_t') >> title('Ajust per minims quadrats a la funcio: C_t=C_0*exp(-K*t)') >> hold off % (b) Vida mitjana. Busquem t=tau t.q. exp(-K*tau)=1/2, i llavors és % clar que: >> tau=log(2)/K tau = 41.7258044011813 % min^(-1) % (c) Càlcul d'una base de polinomis ortogonals fins a grau 2 usant els % temps donats a la taula. % Phi0, Phi1, Phi2: Polinomis de la base canònica de l'ev dels % polinomis de grau <=2: >> Phi0=[0 0 1]% Phi0(t)=1 Phi0 = 0 0 1 >> Phi1=[0 1 0]% Phi1(t)=t Phi1 = 0 1 0 >> Phi2=[1 0 0]% Phi2(t)=t^2 Phi2 = 1 0 0 % Com al problema 6, aquí també farem ortogonalització de Gram- % Schmidt. Per això ens serà útil la funció ddot (definida al % mateix problema com una funció inline) que fa el producte escalar % discret de dos polinomis p i q respecte a les abscisses % x=[x(1),...,x(n)] d'una taula donada. >> ddot=inline('sum(polyval(p,x).*polyval(q,x))','p','q','x'); % Ortogonalització de Gram-Schmidt >> Psi0=Phi0 Psi0 = 0 0 1 >> Psi1=Phi1-ddot(Phi1,Psi0,t)/ddot(Psi0,Psi0,t)*Psi0 Psi1 = 0 1 -75 >> Psi2=Phi2-ddot(Phi2,Psi0,t)/ddot(Psi0,Psi0,t)*Psi0-... > ddot(Phi2,Psi1,t)/ddot(Psi1,Psi1,t)*Psi1 Psi2 = 1 -150 4500 % Llavors, veiem que els polinomis ortogonals buscats de grau <= 2 % són: % Psi0 = [0,0,1] % Psi0(t)=1 % Psi1 = [0,1,-75] % Psi1(t)=t-75 % Psi2 = [1,-150,4500] % Psi2(t)=t^2-150*t+4500 diary off