robert07
Newbie
Offline
Posts: 6
|
 |
«
Embed this message
Reply #1 on: April 13, 2013, 11:56:25 am » |
|
Thanks for asking your frequent question. I am giving you a sample of how to programme an applet of Henon-Heiles? Bellow I the example are a example.Hope that you will be get a benifit from it
Section 1 /* C: osc.c = solve anharmonic oscillator equation */ #include #include #define dt (0.01) #define omega2 (1.0) #define beta (0.4) double f(double x){ return (-omega2*x-beta*pow(x,3)); } void lf2(double x,double p,double *xn,double *pn,double h){ double ph; ph=p+(h/2)*f(x); *xn=x+h*ph; *pn=ph+(h/2)*f(*xn); } void rk4(double x,double p,double *xn,double *pn,double h){ double x1,x2,x3,x4,p1,p2,p3,p4; x1=h*p; p1=h*f(x); x2=h*(p+p1/2); p2=h*f(x+x1/2); x3=h*(p+p2/2); p3=h*f(x+x2/2); x4=h*(p+p3); p4=h*f(x+x3); *xn=x+(x1+2*x2+2*x3+x4)/6; *pn=p+(p1+2*p2+2*p3+p4)/6; } int main(){ double x, p, xn, pn, energy; double d, d1, d2; int solver, n, max=800; FILE *fp; /* 4th order yoshida parameters */ d=pow(2,1.0/3.0); d1=1/(2-d); d2=-d/(2-d); /* choice of scheme */ printf("Which algorithm? (1:Runge-Kutta,2:Leap Frog,3:Yoshida)"); scanf("%d",&solver); /* initial condition */ x=2.0; p=0.0; energy=p*p/2+omega2*x*x/2+beta*pow(x,4)/4; printf("First energy=%10.6f ¥n",energy); /* file open */ switch(solver){ case 1: fp=fopen("rk4.dat","w"); break; case 2: fp=fopen("lf2.dat","w"); break; case 3: fp=fopen("ys4.dat","w"); break; } fprintf(fp,"%10.7f %10.7f¥n",x,p); /* time development */ for(n=1; n switch(solver){ case 1:/* Runge-Kutta */ rk4(x,p,&xn,&pn,dt); x=xn; p=pn; break; case 2:/* Leap Flog */ lf2(x,p,&xn,&pn,dt); x=xn; p=pn; break; case 3:/* Yoshida */ lf2(x,p,&xn,&pn,dt*d1); x=xn; p=pn; lf2(x,p,&xn,&pn,dt*d2); x=xn; p=pn; lf2(x,p,&xn,&pn,dt*d1); x=xn; p=pn; break; } /* write data */ fprintf(fp,"%10.7f %10.7f¥n",x,p); }/* end of time development */ fclose(fp); /* check of the energy conservation */ energy=p*p/2+omega2*x*x/2+beta*pow(x,4)/4; printf("Final energy=%10.6f ¥n",energy); /* end */ return 0; }
Section 2 /* Java: henon.java = Applet program with henon.html:
Henon-Heiles system
code=henon.class width=700 height=600>
*/ import java.applet.Applet; import java.awt.*; import java.awt.event.*; import java.io.*; public class henon extends Applet implements ActionListener, ItemListener, MouseListener{ CheckboxGroup cbg=new CheckboxGroup(); TextField tf=new TextField(15); Button btn1, btn2, btn3; Image img; Graphics bg; String msg; protected int width, height, ON=1, OFF=0; protected int xc=350, yc=350, scale=500; int SW, mx, my; // variables int count, max=1000000; double dt=0.001, t, dd; double x1, x2, p1, p2, energy; double sx, sy; double[][] PS=new double[1000][2];// Poincar¥'e section public void init(){ // canvas width=getSize().width; height=getSize().height; img=this.createImage(width,height); setBackground(Color.white); bg=img.getGraphics(); resetDisplay(); // items Checkbox cb1=new Checkbox("E=1/12",false,cbg); Checkbox cb2=new Checkbox("E=1/8",false,cbg); Checkbox cb3=new Checkbox("E=1/6",false,cbg); add(cb1); cb1.addItemListener(this); add(cb2); cb2.addItemListener(this); add(cb3); cb3.addItemListener(this); add(tf); btn1=new Button("START"); btn2=new Button("STOP"); btn3=new Button("CLEAR"); add(btn1); btn1.addActionListener(this); add(btn2); btn2.addActionListener(this); add(btn3); btn3.addActionListener(this); addMouseListener(this); } public void resetDisplay(){ bg.clearRect(0,60,width,height); bg.setColor(Color.black); bg.drawLine(40,yc,width-40,yc); bg.drawLine(xc,130,xc,height-30); bg.drawImage(img,0,0,this); repaint(); } public void itemStateChanged(ItemEvent ie){ String str=((Checkbox) ie.getItemSelectable()).getLabel(); if(str=="E=1/12"){ energy=1.0/12.0; t=0.0; }else if(str=="E=1/8"){ energy=1.0/8.0; t=0.0; }else if(str=="E=1/6"){ energy=1.0/6.0; t=0.0; } } public void actionPerformed(ActionEvent ae){ if(ae.getSource()==btn1){// START SW=ON; calc(); }else if(ae.getSource()==btn2){// STOP SW=OFF; }else if(ae.getSource()==btn3){// CLEAR SW=OFF; for(int i=0; i<1000; i++){ PS[0]=0.0; PS[1]=0.0; } resetDisplay(); } } public void mousePressed(MouseEvent me){ mx=me.getX( ); my=me.getY( ); if(my<100){ return;// out of bounds }else{ x1=0.0; x2=(double)(mx-xc)/(double)scale; p2=(double)(my-yc)/(double)scale; dd=2*energy-p2*p2-2*V(x1,x2); if(dd>0){ p1=Math.sqrt(dd); tf.setText("x2="+x2+", p2="+p2); SW=ON; }else{ tf.setText("No Good, Try again!"); } } } public void mouseClicked(MouseEvent evt){ } public void mouseReleased(MouseEvent evt){ } public void mouseEntered(MouseEvent evt){ } public void mouseExited(MouseEvent evt){ } public void paint(Graphics g){ int vx, vy; bg.setColor(Color.red); for(int i=0; i<1000; i++){ vx=xc+(int)(scale*PS[0]); vy=yc+(int)(scale*PS[1]); bg.drawRect(vx,vy,1,1); } g.drawImage(img,0,0,this); } // potential energy public double V(double x1, double x2){ return ((x1*x1+x2*x2)/2+x1*x1*x2-x2*x2*x2/3); } // right hand sides public double f1(double x1, double x2, double p1, double p2){ return (p1); } public double f2(double x1, double x2, double p1, double p2){ return (p2); } public double f3(double x1, double x2, double p1, double p2){ return (-x1-2*x1*x2); } public double f4(double x1, double x2, double p1, double p2){ return (-x2-x1*x1+x2*x2); } // for use of Poincar¥'e section public double g1(double x1, double x2, double p1, double p2){ return (1/p1); } public double g2(double x1, double x2, double p1, double p2){ return (p2/p1); } public double g3(double x1, double x2, double p1, double p2){ return ((-x1-2*x1*x2)/p1); } public double g4(double x1, double x2, double p1, double p2){ return ((-x2-x1*x1+x2*x2)/p1); } public void calc(){ double x1_1, x1_2, x1_3, x1_4, x1n; double x2_1, x2_2, x2_3, x2_4, x2n; double p1_1, p1_2, p1_3, p1_4, p1n; double p2_1, p2_2, p2_3, p2_4, p2n; double dx, t_1, t_2, t_3, t_4; if(SW==ON){ count=0; for(int i=0; i x1_1=dt*f1(x1,x2,p1,p2); x2_1=dt*f2(x1,x2,p1,p2); p1_1=dt*f3(x1,x2,p1,p2); p2_1=dt*f4(x1,x2,p1,p2); x1_2=dt*f1(x1+x1_1/2,x2+x2_1/2,p1+p1_1/2,p2+p2_1/2); x2_2=dt*f2(x1+x1_1/2,x2+x2_1/2,p1+p1_1/2,p2+p2_1/2); p1_2=dt*f3(x1+x1_1/2,x2+x2_1/2,p1+p1_1/2,p2+p2_1/2); p2_2=dt*f4(x1+x1_1/2,x2+x2_1/2,p1+p1_1/2,p2+p2_1/2); x1_3=dt*f1(x1+x1_2/2,x2+x2_2/2,p1+p1_2/2,p2+p2_2/2); x2_3=dt*f2(x1+x1_2/2,x2+x2_2/2,p1+p1_2/2,p2+p2_2/2); p1_3=dt*f3(x1+x1_2/2,x2+x2_2/2,p1+p1_2/2,p2+p2_2/2); p2_3=dt*f4(x1+x1_2/2,x2+x2_2/2,p1+p1_2/2,p2+p2_2/2); x1_4=dt*f1(x1+x1_3,x2+x2_3,p1+p1_3,p2+p2_3); x2_4=dt*f2(x1+x1_3,x2+x2_3,p1+p1_3,p2+p2_3); p1_4=dt*f3(x1+x1_3,x2+x2_3,p1+p1_3,p2+p2_3); p2_4=dt*f4(x1+x1_3,x2+x2_3,p1+p1_3,p2+p2_3); x1n=x1+(x1_1+2*x1_2+2*x1_3+x1_4)/6; x2n=x2+(x2_1+2*x2_2+2*x2_3+x2_4)/6; p1n=p1+(p1_1+2*p1_2+2*p1_3+p1_4)/6; p2n=p2+(p2_1+2*p2_2+2*p2_3+p2_4)/6; if((x1<0)&&(x1n>0)){ /* Poincar¥'e section at x1=0 */ dx=-x1; t_1=dx*g1(x1,x2,p1,p2); x2_1=dx*g2(x1,x2,p1,p2); p1_1=dx*g3(x1,x2,p1,p2); p2_1=dx*g4(x1,x2,p1,p2); t_2=dx*g1(x1+dx/2,x2+x2_1/2,p1+p1_1/2,p2+p2_1/2); x2_2=dx*g2(x1+dx/2,x2+x2_1/2,p1+p1_1/2,p2+p2_1/2); p1_2=dx*g3(x1+dx/2,x2+x2_1/2,p1+p1_1/2,p2+p2_1/2); p2_2=dx*g4(x1+dx/2,x2+x2_1/2,p1+p1_1/2,p2+p2_1/2); t_3=dx*g1(x1+dx/2,x2+x2_2/2,p1+p1_2/2,p2+p2_2/2); x2_3=dx*g2(x1+dx/2,x2+x2_2/2,p1+p1_2/2,p2+p2_2/2); p1_3=dx*g3(x1+dx/2,x2+x2_2/2,p1+p1_2/2,p2+p2_2/2); p2_3=dx*g4(x1+dx/2,x2+x2_2/2,p1+p1_2/2,p2+p2_2/2); t_4=dx*g1(x1+dx,x2+x2_3,p1+p1_3,p2+p2_3); x2_4=dx*g2(x1+dx,x2+x2_3,p1+p1_3,p2+p2_3); p1_4=dx*g3(x1+dx,x2+x2_3,p1+p1_3,p2+p2_3); p2_4=dx*g4(x1+dx,x2+x2_3,p1+p1_3,p2+p2_3); t+=(t_1+2*t_2+2*t_3+t_4)/6; x1n=0.0; x2n=x2+(x2_1+2*x2_2+2*x2_3+x2_4)/6; p1n=p1+(p1_1+2*p1_2+2*p1_3+p1_4)/6; p2n=p2+(p2_1+2*p2_2+2*p2_3+p2_4)/6; PS[count][0]=x2n; PS[count][1]=p2n; count+=1; repaint(); }else{ t+=dt; } x1=x1n; x2=x2n; p1=p1n; p2=p2n; } } } }
Section 3 # Perl: kdv.pl = KdV equation according to Zabusky-Kruskal scheme use Tk;
HEIGHT=600;
mw->title("KdV");
mw->Canvas(width=> HEIGHT)->pack(); sub init{ # constants PI=4*atan2(1,1); N=200; N; dt=0.0001; delta=0.022; dt/(3* c_uxxx=( dt/(h**3); # initial condition for(j=0; N+4);j++){ x= j-2); j]=cos( x); j]= j]; } # start canvas->after(10,¥&run); } sub run{ # display canvas->delete('tag'); for( j<=( j++){ j; u[j]*100; canvas->createOval( yy, yy,-tags=>'tag'); } # computations for (1..40) { for( j<=( j++){ c_uux*( j+1]+ j]+ j-1])*( j+1]- j-1]); c_uxxx*( j+2]-2* j+1]+2* j-1]- j-2]); j]= j]- uxxx; } # periodic boundary condition u[ u[1]= N+1]; N+3]= u[ u[4]; unew[ unew[1]= N+1]; N+3]= unew[ unew[4]; # renewal for( j<=( j++){ j]= j]; j]= j]; } } $canvas->after(1, ¥&run); } init(); MainLoop;
Section 4 # Python: BZ.py = 2dim. BZ reaction-diffusion equation import Tkinter import time from numpy import * from math import * # canvas WIDTH=520 HEIGHT=450 canvas=Tkinter.Canvas(width=WIDTH,height=HEIGHT,bg='white') canvas.pack() # model N x M: PBC N, M=30, 60 # parameters pi=4*atan(1.0) dt, h=0.01, 0.1 h2=h*h DX, DY=0.0001, 0.0001 A, B=1.0, 3.0 maxX, maxY=0.0, 0.0 minX, minY=6.0, 6.0 # preparation X=zeros([N+2, M+2],float) Y=zeros([N+2, M+2],float) XX=zeros([N+2, M+2],float) YY=zeros([N+2, M+2],float) rhsX=zeros([N+2, M+2],float) rhsY=zeros([N+2, M+2],float) for j in range(N+2): for k in range(M+2): s, t=pi*j/N, pi*k/M X[j][k]=2.0+0.4*sin(s)*sin(3*t) Y[j][k]=1.0+0.4*sin(s)*sin(t) diagX, offX=1-4*DX*dt/h2, DX*dt/h2 diagY, offY=1-4*DY*dt/h2, DY*dt/h2 # nonlinear reaction terms def SX(u,v): return (A-(B+1)*u+u*u*v) def SY(u,v): return (B*u-u*u*v) # coloring def floatRgb(mag, cmin, cmax): try: # normalize to [0,1] x = float(mag-cmin)/float(cmax-cmin) except: # cmax = cmin x=0.5 blue = min((max((4*(0.75-x), 0.)), 1.)) red = min((max((4*(x-0.25), 0.)), 1.)) green= min((max((4*math.fabs(x-0.5)-1., 0.)), 1.)) return (red, green, blue) def strRgb(mag, cmin, cmax): red, green, blue = floatRgb(mag, cmin, cmax) return "#%02x%02x%02x" % (red*255, green*255, blue*255) # main def run(): global X, XX, Y, YY, rhsX, rhsY, maxX, maxY, minX, minY canvas.delete(Tkinter.ALL) # show color table for j in range (256): color=strRgb(float(j), 0.0, 256.0) canvas.create_rectangle(j+140-6,HEIGHT-40-3,j+140+6,HEIGHT-20+3, width=0,fill=color) # find maxX, maxY, minX, minY for j in range(1, N+1): for k in range(1, M+1): if maxX maxX=X[j][k] if maxY maxY=Y[j][k] if minX>X[j][k]: minX=X[j][k] if minY>Y[j][k]: minY=Y[j][k] # display for j in range(1, N+1): for k in range(1, M+1): x, y=6*(j-1), 20+6*(k-1) # coloring color=strRgb(X[j][k], minX, maxX) canvas.create_rectangle(x+50-3,y-3,x+50+3,y+3,width=0,fill=color) color=strRgb(Y[j][k], minY, maxY) canvas.create_rectangle(x+300-3,y-3,x+300+3,y+3,width=0,fill=color) # run 10 times for i in range(5): for j in range(1, N+1): for k in range(1, M+1): rhsX[j][k]=dt*SX(X[j][k], Y[j][k]) rhsY[j][k]=dt*SY(X[j][k], Y[j][k]) rhsX[j][k]+=offX*(X[j+1][k]+X[j-1][k]+X[j][k-1]+X[j][k+1])+diagX*X[j][k] rhsY[j][k]+=offY*(Y[j+1][k]+Y[j-1][k]+Y[j][k-1]+Y[j][k+1])+diagY*Y[j][k] XX[j][k], YY[j][k]=rhsX[j][k], rhsY[j][k] for j in range(0, N+2): XX[j][0], XX[j][M+1]=XX[j][M], XX[j][1] YY[j][0], YY[j][M+1]=YY[j][M], YY[j][1] for k in range(0, M+2): XX[0][k], XX[N+1][k]=XX[N][k], XX[1][k] YY[0][k], YY[N+1][k]=YY[N][k], YY[1][k] # rename X, Y=XX, YY canvas.after(1, run) run() canvas.mainloop()
Still if you have any question please post here. Thanks
|