problemset3

Z bazant 18366 random walks and diusion problem set 3

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: dt); int blocks=40; //Number of blocks to save out const float pi=3.1415926535897932384626433832795; M. Z. Bazant – 18.366 Random Walks and Diffusion – Problem Set 3 Solutions 13 3500 3000 ρ 2500 2000 1500 Iterations 1000 to 1019 Iterations 1200 to 1219 Iterations 1400 to 1419 1000 500 -4 -3 -2 -1 0 1 2 3 4 x Figure 8: A transient in the density profile, occurring in a much larger version of the simulation. At approximately 1000 iterations a large peak is seen at the top of the standing wave but this dissipates as the simulation progresses. const const const const const const const const int lwalk=int(mx ∗ lrho); //Total initial walkers on LHS int twalk=int(mx ∗ (lrho+rrho)); //Total initial walkers float lspeed=udt∗(1−lrho/rhomax)+ddt; float rspeed=ddt−udt∗(1−rrho/rhomax); float lrate=0.5∗lspeed ∗ lrho; //L. walker intro rate float rrate=0.5∗rspeed ∗ rrho; //R. walker intro rate int mem=twalk∗2; //Total memory float isp=float(blocks)/(2∗mx); //Inverse block spacing inline float rnd() { return (float(rand())+0.5)/RAND MAX; } inline int poisson(float l) { if (l>20) { int r=int(sqrt(−2∗log(rnd())∗l) ∗ cos(2∗pi ∗ rnd())+0.5+l); return r>0?r:0; } else { double a=rnd()∗exp(l)−1,b=1;int r=0; while(a>0) { a−=b ∗ =l/++r; } return r; } } int main() { int a=0,i,j=0,t;float w[mem],v[mem],as[blocks],r,x; while(j++<lwalk) w[a++]=−mx ∗ rnd(); while(j++<twalk) w[a++]=mx ∗ rnd(); M. Z. Bazant – 18.366 Random Walks and Diffusion – Problem Set 3 Solutions for(i=0;i<blocks;i++) as[i]=0; for(t=0;t<iter;t++) { for(i=0;i<a;i++) { j=int((w[i]+mx) ∗ isp); if (j>=0&&j<blocks&&t>100) as[j]++; r=w[i]<dx−mx?lrho ∗ (dx−mx−w[i]):0; r+=w[i]>mx−dx?rrho ∗ (w[i]−mx+dx):0; for(j=0;j<a;j++) { if (abs(w[i]−w[j])<dx) r++; } r/=2∗dx; v[i]=udt∗(1−r/rhomax)+(rand()%2==1?ddt:−ddt); } i=0; while(i<a) { if (abs(w[i]+=v[i])>mx) { v[i]=v[a]; w[i]=w[−−a]; } else i++; } j=poisson(lrate); for(i=0;i<j;i++) w[a++]=lspeed ∗ rnd()−mx; j=poisson(rrate); for(i=0;i<j;i++) w[a++]=mx−rspeed ∗ rnd(); } for(i=0;i<blocks;i++) { r=(float(i)+0.5)/isp−mx; cout << r << &qu...
View Full Document

This note was uploaded on 01/23/2014 for the course MATH 18.366 taught by Professor Martinbazant during the Fall '06 term at MIT.

Ask a homework question - tutors are online