Pendulum Poincare map

liliangbiao

// PendulumPoincare.java
import java.awt.*;
import java.awt.event.*;
import java.awt.Graphics;
public class PendulumPoincare extends Frame
{
  public PendulumPoincare()
  {
  setSize(600,500);
  addWindowListener(new WindowAdapter()
  { public void windowClosing(WindowEvent event)
  { System.exit(0); }}); }
  public void fsystem(double h,double t,double u[],double hf[])
  {
  double a = 0.2, b = 1.0;
  double k = 1.09; double Omega = 0.8;  
  hf[0] = h*u[1];
  hf[1] = h*(-a*u[1]-b*Math.sin(u[0])+k*Math.cos(Omega*t));
  }
  public void map(double u[],int steps,double h,double t,int N)
  {
  double uk[] = new double[N];
  double tk;
  double a[] = { 0.0, 1.0/4.0, 3.0/8.0, 12.0/13.0, 1.0, 1.0/2.0 };
  double c[] = { 16.0/135.0, 0.0, 6656.0/12825.0, 28561.0/56430.0,
                 -9.0/50.0, 2.0/55.0 };
  double b[][] = new double[6][5];
  b[0][0] = b[0][1]= b[0][2] = b[0][3] = b[0][4] = 0.0;
  b[1][0] = 1.0/4.0; b[1][1] = 0.0; b[1][2] = 0.0; b[1][3] = 0.0;
  b[1][4] = 0.0;
  b[2][0] = 3.0/32.0; b[2][1] = 9.0/32.0;
  b[2][2] = 0.0; b[2][3] = 0.0; b[2][4] = 0.0;
  b[3][0] = 1932.0/2197.0; b[3][1] = -7200.0/2197.0;
  b[3][2] = 7296.0/2197.0; b[3][3] = b[3][4] = 0.0;
  b[4][0] = 439.0/216.0; b[4][1] = -8.0;
  b[4][2] = 3680.0/513.0; b[4][3] = -845.0/4104.0; b[4][4] = 0.0;
  b[5][0] = -8.0/27.0; b[5][1] = 2.0; b[5][2] = -3544.0/2565.0;
  b[5][3] = 1859.0/4104.0; b[5][4] = -11.0/40.0;
  double f[][] = new double[6][N];
  for(int i=0;i<steps;i++)
  {
  fsystem(h,t,u,f[0]);
  for(int k=1;k<=5;k++)
  { tk = t+a[k]*h;
  for(int l=0;l<N;l++)
  { uk[l] = u[l];
  for(int j=0;j<=k-1;j++) uk[l] += b[k][j]*f[j][l];
  }
  fsystem(h,tk,uk,f[k]);
  }
  for(int l=0;l<N;l++)
   for(int k=0;k<6;k++) u[l] += c[k]*f[k][l];
  }
  }
  public void paint(Graphics g)
  {
  g.drawRect(40,40,500,450);
  int steps = 1; int N = 2; int counter = 0;
  double h = 0.005;  double Omega = 0.8;  double t = 0.0;
  double u[] = { 0.8, 0.6 };  // initial conditions
  for(int i=0;i<800000;i++)
  {
  t += h; map(u,steps,h,t,N);
  if((counter>10) && (Math.abs(t-2.0*Math.PI*counter/Omega)<0.02))
  {
  int m = ((int)(10.0*u[0])+350);  int n = ((int)(40.0*u[1])+200);
  g.drawLine(m,n,m,n);
  }
  if(t > 2.0*Math.PI*counter/Omega) { counter++; }
  }
  }
  public static void main(String[] args)
  {
  Frame f = new PendulumPoincare(); f.setVisible(true);
  }
}

hsfy919

钟摆系统的庞加莱截面程序，感谢liliangbiao分享，只是由于语言的差异看不太懂，希望有人出来解释一下，不知道与论坛上已有的matlab程序在原理上有什么不同？

雨人

liliangbiao先生用java写的，map函数不太明白，太多赋值语句了，如果做映射应该在paint函数里面吧，但是不太清楚，还是请别的高手路过吧！

喳喳123

没有学过JAVA呢，看您的程序，只能了解您编程的思想，但是还是很受用的。试着用MATLAB实现您这一思想。

liliangbiao

这就对了，只需要了解编程思想就可以了，试试用Matlab或者其他语言编写一下，应该不难理解的！编写对了，欢迎共享！

lalama

调试的时候，我只有半个图，是不是我开的窗口太小了？