Rootspoly24

Z Skrypty dla studentów Ekonofizyki UPGOW

Wersja Marcin (dyskusja | edycje) z dnia 12:36, 17 lut 2011

(różn.) ← poprzednia wersja | przejdź do aktualnej wersji (różn.) | następna wersja → (różn.)

Zmiana bazy i iloczyn skalarny

Mamy \((x,y)=\sum_{i=1}^{N}x_i y_j\)

b1=[1;1;1]
b2=[1;0;1]
b3=[1;0;-1]
C=[b1,b2,b3]
C'*C
M=[b1'*b1,b1'*b2,b1'*b3;b2'*b1,b2'*b2,b2'*b3;b3'*b1,b3'*b2,b3'*b3]

Ortogonalizacja Grama-Schmidta

e1=b1
e2=b2-(e1'*b2)/(e1'*e1)*e1
e3=b3-( (e1'*b3)/(e1'*e1)*e1 + (e2'*b3)/(e2'*e2)*e2 )

E=[e1,e2,e3]
quiver3([0,0,0],[0,0,0],[0,0,0], C(1,:), C(2,:), C(3,:), 0,'r')
hold on
quiver3([0,0,0],[0,0,0],[0,0,0], E(1,:), E(2,:), E(3,:), 0,'b')

Ciekawa krzywa

\(x^2+(y-\sqrt{x^2})^2=1\)

Zera wielomianów 24-tego stopnia

Miejsca zerowe na płaszczyźnie zespolonej, wszystkich wielomianiów stopnia 24 ze współczynnikami +1 lub -1.

Mamy \(24*2^24\) punków. Można je zhistogramować.

Roots

close all; 
clear all;
 
subplot(1,2,1)
xlim([-2,2])
ylim([-2,2])
hold on; 
 
subplot(1,2,2)
xlim([-2,2])
ylim([-2,2])
hold on; 
 
tic; 
 
 
 
for i=1:11100; 
  x=roots( (2*(rand(1,24)>0.5))-1 );
  x2=roots( 2*rand(1,24)-1 );
  subplot(1,2,1)
  plot(real(x),imag(x),'.');
  subplot(1,2,2)
  plot(real(x2),imag(x2),'r.');
 # title(num2str(i))
  if (mod(i,1000)==0) 
    drawnow(); 
  end
end;
 
toc
imagesc(-2:2,-2:2,log(a+1)')

close all; 
clear all;
 
 
nx=1024;
ny=1024; 
N=[nx,ny];
Min=[-2,-2];
Max=[ 2, 2];
 
a=zeros(nx,ny);
tic; 
for i=1:122000; 
  x=roots( (2*(rand(1,24)>0.5))-1 );
  ix=floor ( (( real(x)-Min(1))./(Max(1)-Min(1))).*(N(1)-1)) + 1;
  iy=floor ( (( imag(x)-Min(2))./(Max(2)-Min(2))).*(N(2)-1)) + 1;
  a(ix+(iy-1)*nx)+=1;
end;
toc

close all; 
clear all;
n=1000000
 
nth=24;
nx=1024;
ny=1024; 
N=[nx,ny];
Min=[-2,-2];
Max=[ 2, 2];
 
a=zeros(nx,ny);
tic; 
for i=1:n; 
  x=roots( (2*(rand(1,nth)>0.5))-1 );
  ix=floor ( (( real(x)-Min(1))./(Max(1)-Min(1))).*(N(1)-1)) + 1;
  iy=floor ( (( imag(x)-Min(2))./(Max(2)-Min(2))).*(N(2)-1)) + 1;
  a(ix+(iy-1)*nx)+=1;
  if (mod(i,5000)==0) 
    norma=sum(sum(a))/(nx*ny);
    image(-2:2,-2:2,(10/norma*log(a+1))')
    drawnow(); 
  endif
end;
toc/n*2^nth/3600
norma=sum(sum(a))/(nx*ny);
image(-2:2,-2:2,(10/norma*log(a+1))')

x=linspace(-10,10,200);
data=normrnd (2.1,pi,350,1);
hist (  data , -10:0.5:10,2 ) 
hold on 
x=linspace(-15,15,200); 
plot(x,normpdf(x,2.1,pi),'g')
plot(x,normpdf(x,mean(data),sqrt(cov(data))),"r-")
plot (x',mvnpdf (x',mean(data),cov(data)), "+")