Rootspoly24

Z Skrypty dla studentów Ekonofizyki UPGOW

(Różnice między wersjami)

Skocz do: nawigacji, wyszukiwania

Aktualna wersja na dzień 16:35, 24 lut 2011

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

#p=[1,1.2,.2,1,4]
clear all
p=[1 rand(1,14)];
N=length(p)-1;
 
x0= 0.1*( (rand(N,1)-0.5)+I*(rand(N,1)-0.5))   ;
 
x0=roots(p)+rand(N,1)*0.001;
 
  plot(x0,'*r');
  drawnow();
  hold on
 
 
for i = 1:11
 
  for i = 1:N;
    pr(i)= prod( circshift(x0,1-i)(1)-circshift(x0,1-i)(2:N) ); 
  endfor
  x0=x0-polyval(p,x0)./pr';
 
  plot(x0,'.');
  drawnow();
 
endfor
 
hold off

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)), "+")