Physical Simulation

Physical Simulation Example

The following example shows how you might integrate a numerical simulation with a Java® animation. Here, we have considered a simple spring-mass-damper system of the form

    d2x/dt2 + 2xw0dx/dt + w02 x = 0

and computed the solution numerically using an explicit Forward Euler time-stepping scheme. Try increasing the size of the time step and notice that the simulation becomes unstable when the growth factors become larger than 1. When using explicit time integration, we must therefore choose our time step to be smaller than the critical time step.

Now try using an implicit Backward Euler time-stepping scheme. In this case, the simulation remains stable even for large time steps because the growth factors are always smaller than 1. You may also wish to determine the exact solution to the differential equation and compare it to the numerical solution.

import java.awt.*;
import java.awt.event.*;
import javax.swing.*;

class Vec2D {
    private float[] vec;

    public Vec2D(float fX, float fV) {
        vec = new float[2];
        vec[0] = fX;
        vec[1] = fV;
    }

    public void translate(float fDx) {
        vec[0] += fDx;
    }

    public void setPos(float fX) {
        vec[0] = fX;
    }

    public void setVel(float fV) {
        vec[1] = fV;
    }

    public float getPos() {
        return vec[0];
    }

    public float getVel() {
        return vec[1];
    }
}

class Matrix2D {
    private float[][] mat;

    public Matrix2D(float a11, float a12, float a21, float a22) {
        mat = new float[2][2];
        mat[0][0] = a11;
        mat[0][1] = a12;
        mat[1][0] = a21;
        mat[1][1] = a22;
    }

    public void multiply(Vec2D vec) {
        float fX = mat[0][0] * vec.getPos() + mat[0][1] * vec.getVel();
        float fV = mat[1][0] * vec.getPos() + mat[1][1] * vec.getVel();
       vec.setPos(fX);
       vec.setVel(fV);
    }

    public void invert() {
        float det = mat[0][0]*mat[1][1]-mat[0][1]*mat[1][0];
        float tmp = mat[0][0];
       mat[0][0] = mat[1][1] / det;
       mat[1][1] = tmp / det;
       mat[0][1] = -mat[0][1] / det;
       mat[1][0] = -mat[1][0] / det;
    }
}

class Ball {
    Vec2D mXState;  // x position and velocity.
    Vec2D mYState;  // y position and velocity.
    Matrix2D mMatrix;
    int miRadius;
    int miWindowWidth, miWindowHeight;

    public Ball(int iRadius, int iW, int iH) {
        miRadius = iRadius;
        mXState = new Vec2D(0.0f, 0.0f);
        mYState = new Vec2D(0.0f, 0.0f);
        miWindowWidth = iW;
        miWindowHeight = iH;
    }

    public void setPosition(float fXPos, float fYPos) {
        mXState.setPos(fXPos);
        mYState.setPos(fYPos);
    }

    public void setVelocity(float fXVel, float fYVel) {
        mXState.setVel(fXVel);
        mYState.setVel(fYVel);
    }

    public void setParams(float fXi, float fW0, float fDt, boolean explicit) {
        float fReal1 = 0.0f, fImag1 = 0.0f;  // First eigenvalue.
        float fReal2 = 0.0f, fImag2 = 0.0f;  // Second eigenvalue.
        float G1, G2;  // Growth factors.

        // Determine the eigenvalues.
        if (fXi < 1.0f) {
            fReal1 = fReal2 = -fW0*fXi;
            fImag1 = (float)(fW0*Math.sqrt(1-fXi*fXi));
            fImag2 = -fImag1;
            System.out.println("System is underdamped.");
            System.out.println("Eigenvalues are: " + fReal1 + " +/- " + fImag1 + "i");
        }
        else {
            fReal1 = -fW0*(fXi + (float)Math.sqrt(fXi*fXi-1));
            fReal2 = -fW0*(fXi - (float)Math.sqrt(fXi*fXi-1));
            System.out.println("System is overdamped or critically damped.");
            System.out.println("Eigenvalues are: " + fReal1 + " and " + fReal2 + "i");
        }

        if (explicit) {
            // Forward Euler.
            mMatrix = new Matrix2D(1.0f, fDt, -fW0*fW0*fDt, 1-2*fXi*fW0*fDt);

           G1 = (float)Math.sqrt(Math.pow(1+fReal1*fDt,2.0) + Math.pow(fImag1*fDt,2.0));
           G2 = (float)Math.sqrt(Math.pow(1+fReal2*fDt,2.0) + Math.pow(fImag2*fDt,2.0));
        }
        else {
            // Backward Euler.
            mMatrix = new Matrix2D(1.0f, -fDt, fW0*fW0*fDt, 1+2*fXi*fW0*fDt);
            mMatrix.invert();

            G1 = (float)(1.0/Math.sqrt(Math.pow(1-fReal1*fDt,2.0) + Math.pow(fImag1*fDt,2.0)));
            G2 = (float)(1.0/Math.sqrt(Math.pow(1-fReal2*fDt,2.0) + Math.pow(fImag2*fDt,2.0)));
        }
        System.out.println("Growth factors are " + G1 + " and " + G2);
    }

    public int getXPos() {
        return (int)mXState.getPos();
    }

    public int getYPos() {
        return (int)mYState.getPos();
    }

    public void draw(Graphics g) {
        g.setColor(Color.red);
        g.fillOval(miWindowWidth/2+(int)mXState.getPos()-miRadius,
                        miWindowHeight/2+(int)mYState.getPos()-miRadius,
                        2*miRadius, 2*miRadius);
    }

    // Update the position of the ball.
    void move() {
        mMatrix.multiply(mYState);
    }
}

public class Animation extends JApplet implements Runnable, ActionListener {
    int miFrameNumber = 0;
    int miTimeStep;
    Thread mAnimationThread;
    boolean mbIsPaused = false;
    Button mButton;
    Ball ball;

    public void init() {
        // Time step in milliseconds.
        miTimeStep = 20; // Try changing this to (a) 50 ms and (b) 60 ms.

        // Initialize the parameters of the ball.  The parameters refer to the
        // differential equation:  d^2 x/dt^2 + 2 xi w0 dx/dt + w0^2 x = 0

        int iRadius = 15;
        float fXPos = 0.0f;      // Initial x displacement
        float fYPos = 100.0f;    // Initial y displacement
        float fXVel = 0.0f;      // Initial x velocity
        float fYVel = 0.0f;      // Initial y velocity
        float fXi = 0.05f;       // xi
        float fW0 = 2.0f;        // w0
        boolean explicit = true; // true: forward Euler, false: backward Euler

        ball = new Ball(iRadius, getSize().width, getSize().height);
        ball.setPosition(fXPos, fYPos);
        ball.setVelocity(fXVel, fYVel);
        ball.setParams(fXi, fW0, miTimeStep/1000.0f, explicit);

        // Create a button to start and stop the animation.
        mButton = new Button("Stop");
        getContentPane().add(mButton, "North");
        mButton.addActionListener(this);

        // Create a JPanel subclass and add it to the JApplet.  All drawing
        // will be done here, do we must write the paintComponent() method.
        // Note that the anonymous class has access to the private data of
        // class Animation, because it is defined locally.
        getContentPane().add(new JPanel() {
            public void paintComponent(Graphics g) {
                // Paint the background.
                super.paintComponent(g);

                // Display the frame number.
                g.drawString("Frame " + miFrameNumber, getSize().width/2 - 40,
                                      getSize().height - 15);

                // Draw the rubber band.
                g.drawLine(getSize().width/2, 0,
                getSize().width/2+ball.getXPos(),
                getSize().height/2+ball.getYPos());

                // Draw the ball.
                ball.draw(g);
            }
        }, "Center");
    }

    public void start() {
        if (mbIsPaused) {
            // Don't do anything.  The animation has been paused.
        } else {
            // Start animating.
            if (mAnimationThread == null) {
                mAnimationThread = new Thread(this);
            }
            mAnimationThread.start();
        }
    }

    public void stop() {
        // Stop the animating thread by setting the mAnimationThread variable
        // to null.  This will cause the thread to break out of the while loop,
        // so that the run() method terminates naturally.
        mAnimationThread = null;
    }

    public void actionPerformed(ActionEvent e) {
        if (mbIsPaused) {
            mbIsPaused = false;
            mButton.setLabel("Stop");
            start();
        } else {
            mbIsPaused = true;
            mButton.setLabel("Start");
            stop();
        }
    }

    public void run() {
        // Just to be nice, lower this thread's priority so it can't
        // interfere with other processing going on.
        Thread.currentThread().setPriority(Thread.MIN_PRIORITY);

        // Remember the starting time.
        long startTime = System.currentTimeMillis();

        // Remember which thread we are.
        Thread currentThread = Thread.currentThread();

        // This is the animation loop.
        while (currentThread == mAnimationThread) {
            // Draw the next frame.
            repaint();

            // Advance the animation frame.
            miFrameNumber++;

            // Update the position of the ball.
            ball.move();

            // Delay depending on how far we are behind.
            try {
                startTime += miTimeStep;
                Thread.sleep(Math.max(0,
                             startTime-System.currentTimeMillis()));
            }
            catch (InterruptedException e) {
                break;
            }
        }
    }
}