stick2

int fRate = 30;
float dt = 1.0/(float)fRate;
 
float E = 1.0; // Young's modulus [Pa]
float r = 0.01; // radius
float K = r*0.5; // moment of inertia
float rho = 0.1; // density[kg/m^3]
final int N = 200; // number of section
float dx = 0.1; // length of section
 
float[] y = new float[N];
float[] y_temp = new float[N];
float[] dy = new float[N];
 
void setup(){
  frameRate(fRate);
  size(500,400);
  noFill();
  stroke(100,50,20);
  strokeWeight(8);
 
  //initial condition
  for(int i=0;i<N;i++){
    y[i] = 0;
    dy[i] = 0;
    y_temp[i] = 0;
  }
}
 
void draw(){
  int n=0;
  while(n++<200){
    calc();
  }
 
  background(255);
 
  beginShape();
  for(int i=0;i<N;i++){
    curveVertex(i*2 + 20,y[i] + 200);
  }
  endShape();
}
 
// calculation of powers
void calc(){
  for(int i=2;i<N-2;i++){
    dy[i] = dy[i] - (y[i-2] -4*y[i-1] + 6*y[i] - 4*y[i+1] + y[i+2])
      *E*K*K*dt*dt/(rho*dx*dx*dx*dx);
    y_temp[i] = y[i] + dy[i];
  }
 
  for(int i=0;i<N;i++){
    y[i] = y_temp[i];
  }
 
  y[0] = 0;
  y[N-1] = 0;
}
 
// mousePressed Event
void mousePressed(){
  if(mouseX>=20 && mouseX<= 2*N + 20){
    int n = round((mouseX-20)/2);
    dy[n] = 0.1;
  }
}