fork download
copy 
  1.  
  2. #include <stdlib.h>
  3. #include <conio.h>
  4. #include <math.h>
  5. #include <iostream>
  6. #include <fstream>
  7. #include <iomanip>
  8. #include <cmath>
  9. #include <string>
  10. #include <complex>
  11. #include <vector>
  12. #define PI 3.141592654
  13.  
  14. using namespace std;
  15.  
  16. /* function prototypes */
  17. complex<double> dnx(double, vector<complex<double> > , vector<complex<double> >, int );
  18. complex<double> rk4_dn1(complex<double>(*)(double, vector<complex<double> >, vector<complex<double> >, int ),
  19.                double, double, vector<complex<double> >, vector<complex<double> >, int);
  20.                double sampleNormal ();
  21.  
  22.  
  23. /* global variables */
  24. const double g = 0;               
  25. const double r = -1;
  26. const double u = 0.00;                
  27. const double w = 100.0;
  28.  
  29. const complex<double> i (0,1);
  30. int main()
  31.  
  32.  
  33. {
  34.     int  n;
  35.  
  36.     int tmax;
  37. 	int tt = 5000;                   // number of first-order equations 
  38.     double ti, tf, dt, sigma,mu,z,q,N ;
  39.     vector<complex<double> > xi, xf;
  40.     double  j;
  41.     int i, y,d;
  42.     int m=0;
  43.  
  44.  
  45. /* output: file and formats */
  46.  
  47.  
  48.     ofstream file;
  49.     file.open ("provavet11(om100(2g))).dat"); 
  50. 	       // write results to this file
  51. /* output format */
  52.     file.precision(6);
  53.     file.setf(ios::scientific | ios::showpoint);
  54.     cout.precision(6);
  55.     cout.setf(ios::fixed | ios::showpoint);
  56.  
  57.  
  58.  
  59.  
  60.  
  61.     printf("Number of equations\n");
  62.     scanf("%d", &n);
  63.     printf("Particles in central cav.\n");
  64.     scanf("%d", &N);
  65.     printf("Sigma\n");
  66.     scanf("%d", &q);
  67.  
  68. /* initial information */
  69.  
  70.     xi.reserve(n);
  71.     xf.reserve(n);
  72.  
  73.  
  74. 	    // initial value for variable t
  75.     for(y=0; y<=n-1; y++)
  76.  
  77. 	{
  78.  
  79.     //scanf("%f\n", xi[y]);
  80.     //cin >> xi[y];
  81.     //   }
  82.         // step size for integration (s)
  83.     xi[y]= N* exp(-pow((y-(n - 1) /2.),2)/(2 * q * q));
  84. 	cout << xi[y] << "\n";
  85.  
  86.  
  87.          }
  88.  
  89. 			   ti = 0.0;
  90. 			   dt = 0.00005;        			  
  91. 			   tmax = 4000;          // integrate till tmax (s)
  92.  
  93. /* Gaussian randomly distributed array 
  94. for (tt = 0; i <= 10000; tt+dt) {
  95.  
  96. 	eta[tt]= sampleNormal();
  97.  
  98. }
  99.  
  100.  /*
  101.  
  102. /* integration of ODE */
  103.  
  104.  
  105.  
  106.     while (ti <= tmax)
  107.     {
  108.  
  109.     	/*fracpart= modf(ti, &j);*/
  110.  
  111. 		m=m+1;
  112.  
  113. 		y==0;
  114.  
  115.          /* z = real(xi[2*y]*xi[2*y+1]);*/
  116.  
  117.          // fprintf(fp, "%f", xi[y]);
  118.  
  119.           d=(n-2)/2;
  120. 		  if( m%20000 == 0)
  121.  
  122. 		  {
  123.  
  124.           for (y = 0; y <= d; y++) {
  125.  
  126.               file <<setw(8)<< real(xi[2*y]*xi[2*y+1])+imag(xi[2*y]*xi[2*y+1])<< setw(8)<< "";
  127.  
  128.  
  129.             } 
  130.             file << endl;
  131.        }
  132.  
  133.  
  134.  
  135.  
  136.      tf = ti + dt;
  137. /*=================================*/
  138.      rk4_dn1(dnx, ti, tf, xi, xf, n);
  139. /*=================================*/
  140. /* prepare for the next step*/
  141.         ti = tf;
  142.         for (i = 0; i<=n-1; i = i+1)
  143.         { 
  144.            xi[i] = xf[i];
  145.         } 
  146.  
  147.  
  148.  
  149.    }
  150.     system ("pause");
  151.     return 0;
  152. }
  153.  
  154. /*============================================================== 
  155.   System of first order differential equations for the RK solver
  156.   
  157.   For a system of n first-order ODEs
  158.   x [] array - x values
  159.   dx[] array - dx/dt values
  160.   
  161.   
  162. ==============================================================*/
  163.  
  164.     complex<double> dnx(double t, vector<complex<double> > x, vector<complex<double> > dx, int n)
  165. {
  166. 		const complex<double> i (0,1);
  167. 		int j,d;
  168. 		d=(n-2)/2;
  169.  
  170. /* first cav */
  171.       dx[0] = -i*w*x[0]-i*g*abs(x[0])*x[1]-i*r*(x[n-2]+x[2]-2.0*x[0])-(1.0-i)*u*sampleNormal()*x[0];
  172.       dx[1] = i*w*x[1]+i*g*abs(x[1])*x[0]+i*r*(x[n-1]+x[3]-2.0*x[1])-(1.0+i)*u*sampleNormal()*x[1];
  173.       dx[n-2] = -i*w*x[n-2]-i*g*abs(x[n-2])*x[n-1]-i*r*(x[n-4]+x[0]-2.0*x[n-2])-(1.0-i)*u*sampleNormal()*x[n-2];
  174.       dx[n-1] = i*w*x[n-1]+i*g*abs(x[n-1])*x[n-2]+i*r*(x[n-3]+x[1]-2.0*x[n-1])-(1.0+i)*u*sampleNormal()*x[n-1];
  175.  
  176.  
  177. for(j=1;j<=d-1;j++){
  178.  
  179.  
  180.  
  181.       dx[2*j] = -i*w*x[2*j]-i*g*abs(x[2*j])*x[2*j+1]-i*r*(x[2*j-2]+x[2*j+2]-2.0*x[2*j])-(1.0-i)*u*sampleNormal()*x[2*j];
  182.       dx[2*j+1] = i*w*x[2*j+1]+i*g*abs(x[2*j+1])*x[2*j]+i*r*(x[2*j-1]+x[2*j+3]-2.0*x[2*j+1])-(1.0+i)*u*sampleNormal()*x[2*j+1];
  183. }
  184.  
  185.  
  186. }    
  187.  
  188. /*==========================================================
  189.       Runge-Kutta 4th-order
  190. ----------------------------------------------------------
  191.  call ...
  192.  dnx(t,x[],dx[],n)- functions dx/dt   
  193.  input ...
  194.  ti    - initial time
  195.  tf    - solution time
  196.  xi[]  - initial values 
  197.  n     - number of first order equations
  198.  output ...
  199.  xf[]  - solutions
  200. ==========================================================*/
  201. complex<double> rk4_dn1(complex<double>(dnx)(double, vector<complex<double> >, vector<complex<double> >, int), 
  202.                double ti, double tf, vector<complex<double> > xi, vector<complex<double> > xf, int n)
  203. {
  204.      vector<complex<double> >  x, dx;
  205.       vector<complex<double> > k1,k2,k3,k4;
  206.       double t,h;
  207.  
  208.       int j;
  209.  
  210.       h = tf-ti;
  211.       t = ti;
  212. //k1
  213.       dnx(t, xi, dx, n);
  214.       for (j = 0; j<=n-1; j = j+1)
  215.         {
  216.           k1[j] = h*dx[j];
  217.           x[j]  = xi[j] + k1[j]/2.0;  
  218.         }      
  219. //k2
  220.  
  221.       dnx(t+h/2.0, x, dx, n);
  222.       for (j = 0; j<=n-1; j = j+1)
  223.         {
  224.           k2[j] = h*dx[j];
  225.           x[j]  = xi[j] + k2[j]/2.0;  
  226.         }
  227. //k3
  228.       dnx(t+h/2.0, x, dx, n);
  229.       for (j = 0; j<=n-1; j = j+1)
  230.         {
  231.           k3[j] = h*dx[j];
  232.           x[j]  = xi[j] + k3[j];  
  233.         }      
  234. //k4 and result      
  235.       dnx(t+h, x, dx, n);
  236.       for (j = 0; j<=n-1; j = j+1)
  237.         {
  238.           k4[j] = h*dx[j];
  239.           xf[j] = xi[j] + k1[j]/6.0+k2[j]/3.0+k3[j]/3.0+k4[j]/6.0;
  240.         }      
  241.  
  242.  
  243.  
  244.  
  245.     return 0.0;
  246.  
  247. }
  248. /* random gaussian generator */
  249. double sampleNormal() 
  250. {
  251.     double u = ((double) rand() / (RAND_MAX)) * 2 - 1;
  252.     double v = ((double) rand() / (RAND_MAX)) * 2 - 1;
  253.     double r = u * u + v * v;
  254.     if (r == 0 || r > 1) return sampleNormal();
  255.     double c = sqrt(-2 * log(r) / r);
  256.     return u * c;
  257. }
  258.  
  259.  
Compilation error #stdin compilation error #stdout 0s 0KB
comments (?)
stdin
copy 
Standard input is empty
compilation info 
prog.cpp:3:19: fatal error: conio.h: No such file or directory
 #include <conio.h>
                   ^
compilation terminated.
stdout
copy 
Standard output is empty
https://ideone.com/Eckq3Z
language:
C++ (gcc 8.3)
created:
9 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!