Metoda Monte Carlo jest losowa metoda numeryczna. Program znajduje w sposob losowy najlepsze parametry planet. Metoda ta nie jest zbyt wydajna. Polega na tym, ze losujemy zespol parametrow i sprawdzamy jaka wartosc ma funkcja dopasowania, po czym znowu losujemy parametry i sprawdzamy. Jezeli dopasowanie jest lepsze zachowujemy je w pamieci. Czynimy tak do czasu znalezienia odpowiednio dobrego dopasowania lub po wykonaniu odpowiedniej liczby krokow. W programie nalezy podac zakrs parametrow (granice strzalu) oraz nazwe pliku wyjsciowego. -------------------- <verbatim> /*MonteCarlo*/ #include<stdio.h> #include<math.h> #include<stdlib.h> #define LKMAX 100000000 #define LP 11 #define M 142 #define N 2 #define KMAX 150 #define DZIEN 86400 #define PMAX 1000 #define VMAX 100 #define PI 3.1415926 #define B 0.0001 float kepler(float okres,float tau, float t, float e, float b); main() { FILE *dane; FILE *wynik; char *plik_dane; char *plik_wynik; float t[M],v[M],d[M]; int i,j,lk; float pom,v0,chi,chi2,sum,EE,niu; float K[N],P[N],e[N],w[N],tau[N]; int pomocnicza; plik_dane="hr7.dat"; dane=fopen(plik_dane,"r"); for (j=0; j<=M-1; ++j) { fscanf(dane,"%f%f%f\n",&t[j],&v[j],&d[j]); t[j]=t[j]*DZIEN; } fclose(dane); /* losowanie parametrow poczatkowych */ pom=rand(); v0=VMAX*(pom/RAND_MAX-0.5); for (i=0; i<N; ++i) { pom=rand(); K[i]=KMAX*pom/RAND_MAX; pom=rand(); e[i]=pom/RAND_MAX; pom=rand(); w[i]=2*PI*pom/RAND_MAX; pom=rand(); P[i]=PMAX*DZIEN*pom/RAND_MAX; pom=rand(); tau[i]=P[i]*pom/RAND_MAX+t[0]; } /* obliczamy chi^2/(M-LP-1) */ chi=0.0; for (j=0; j<=M-1; ++j) { sum=0.0; for (i=0; i<N; ++i) { EE=kepler(P[i],tau[i],t[j],e[i],B); niu=2*atan(sqrt((1+e[i])/(1-e[i]))*tan(EE/2));; if (niu<0) niu=niu+2*PI; sum=sum+K[i]*(cos(w[i]+niu)+e[i]*cos(w[i])); } chi=chi+((v[j]-v0-sum)/d[j])*((v[j]-v0-sum)/d[j]); } chi=chi/(M-LP-1); /* zapisujemy parametry oraz chi do pliku */ plik_wynik="fit9.dat"; wynik=fopen(plik_wynik,"a"); for (i=0; i<N; ++i) { fprintf(wynik,"%10.4f%10.4f%10.4f%15.4f%15.4f",K[i],e[i],w[i]*180/PI,P[i]/DZIEN,tau[i]/DZIEN); } fprintf(wynik,"%10.4f%10.4f\n",v0,chi); /* losujemy parametry i sprawdzamy, czy daja lepsze dopasowanie. Jesli tak, to zapisaujemy dane do pliku */ pomocnicza=0; for (lk=1; lk<=LKMAX; ++lk) { pom=rand(); v0=VMAX*(pom/RAND_MAX-0.5); for (i=0; i<N; ++i) { pom=rand(); K[i]=KMAX*pom/RAND_MAX; pom=rand(); e[i]=pom/RAND_MAX; pom=rand(); w[i]=2*PI*pom/RAND_MAX; pom=rand(); P[i]=PMAX*DZIEN*pom/RAND_MAX; pom=rand(); tau[i]=P[i]*pom/RAND_MAX+t[0]; } chi2=0.0; for (j=0; j<=M-1; ++j) { sum=0.0; for (i=0; i<N; ++i) { EE=kepler(P[i],tau[i],t[j],e[i],B); niu=2*atan(EE/2); if (niu<0) niu=niu+2*PI; sum=sum+K[i]*(cos(w[i]+niu)+e[i]*cos(w[i])); } chi2=chi2+((v[j]-v0-sum)/d[j])*((v[j]-v0-sum)/d[j]); } chi2=chi2/(M-LP-1); if (chi2<chi) { chi=chi2; for (i=0; i<N; ++i) { fprintf(wynik,"%10.4f%10.4f%10.4f%15.4f%15.4f",K[i],e[i],w[i]*180/PI,P[i]/DZIEN,tau[i]/DZIEN); } fprintf(wynik,"%10.4f%10.4f\n",v0,chi); } pomocnicza=pomocnicza+1; if (pomocnicza==100000) { fclose(wynik); wynik=fopen(plik_wynik,"a"); pomocnicza=0; } } fclose(wynik); return 0; } /* funkcja rozwiazujaca rownanie Keplera */ float kepler(float okres,float tau, float t, float e, float b) { float MM,EE,f0,f1,f2,f3,d1,d2,d3,pi; pi=3.1415926; MM=2*pi*(t-tau)/okres; EE=MM; f0=EE-e*sin(EE)-MM; f1=1-e*cos(EE); f2=e*sin(EE); f3=e*cos(EE); d1=-f0/f1; d2=-f0/(f1+d1*f2/2); d3=-f0/(f1+d2*f2/2+d2*d2*f3/6); while ((d3>b) || (d3<-b)) { EE=EE+d3; f0=EE-e*sin(EE)-MM; f1=1-e*cos(EE); f2=e*sin(EE); f3=e*cos(EE); d1=-f0/f1; d2=-f0/(f1+d1*f2/2); d3=-f0/(f1+d2*f2/2+d2*d2*f3/6); } return EE; } </verbatim> -- Main.CezaryMigaszewski - 26 Feb 2005 ---- Brakuje plik wejść _hr7.dat_ - Też mogłybyś zrobić test żeby uważać użytkownik jeśli plik nie istnieje - jest to więcej ,,user-friendly" niż błęd _Naruszenie ochrony pamięci_ ;) np: <verbatim> plik_dane="hr7.dat"; dane=fopen(plik_dane,"r"); if(dane==NULL) { printf("Niestety, nie widze plik hr7.dat . Prosze go przygotowac.\n"); exit (0); }; </verbatim>
This topic: Main
>
TWikiUsers
>
CezaryMigaszewski
>
MonteCarlo
Topic revision:
15 Mar 2005,
BoudRoukema
(raw view)
Copyright © CC-BY-SA by the contributing authors. All material on this collaboration platform is copyrighted under CC-BY-SA by the contributing authors unless otherwise noted.
Ideas, requests, problems regarding Foswiki?
Send feedback