Runge-Kutta method for Simple pendulum in java

[sefiroths]

Hi, I'm trying to study Runge-Kutta method and apply on a simple pendulum.
Using a timestep dt=0.1 (h=0.1) the pendulum increses energy of 1% every period...
while in this site:
https://www.myphysicslab.com/pendulum/pendulum/pendulum-en.html
decreses energy about 0.01%

What am I doing wrong?
There is a software where can I compare my camculated values?
thanks

import java.awt.BasicStroke;
import java.awt.Color;
import java.awt.Graphics;
import java.awt.Graphics2D;
import java.awt.event.MouseEvent;
import java.awt.event.MouseMotionListener;
import java.awt.geom.Point2D;
import java.util.ArrayList;

import javax.swing.JFrame;
import javax.swing.JPanel;class Sistema{
    double g=10;//9.81;
    double m1 = 1;
    double th1 = Math.PI/2;
    double w1=0;
    double L1=1;

    double dt=0.001;

    final int N = 2;

    double energiaIniziale, yPotenziale0;

    ArrayList<Point2D.Double> tracciato = new ArrayList<Point2D.Double>();

    public Sistema(){
        double y0 = getY();

        yPotenziale0 = m1*L1;
        energiaIniziale = energiaCinetica()+energiaPotenziale();

        System.out.println("y0:"+y0);
        System.out.println("y:"+yPotenziale0 + " Tot."+energiaIniziale+" = C:"+energiaCinetica()+" P:"+energiaPotenziale());

    }

    synchronized ArrayList<Point2D.Double> getTracciato(){
        return tracciato;
    }

    synchronized double energiaCinetica(){
        return m1*L1*L1*w1*w1/2;
    }

    synchronized double energiaPotenziale(){
        return m1*g*yPotenziale0+g*m1*getY();
    }

    synchronized double getX(){
        double x1 = L1*Math.sin(th1);
        return x1;
    }

    synchronized double getY(){
        double y1 = -L1*Math.cos(th1);
        return y1;
    }
    void simulate(){
      
        synchronized(this){
            double yin[] = new double[N];
            yin[0] = th1;
            yin[1] = w1;
            double yout[] = runge_kutta(0, yin, dt);
            this.th1 = yout[0];
            this.w1 = yout[1];
          
            //System.out.println(""+Math.toDegrees(th1)+" "+w1);

            tracciato.add(new Point2D.Double(getX(), getY()));
        }
    }

    double[] derivs(double xin, double yin[])
    {
        /*
            yin[0] = th1[i];
            yin[1] = w1[i];
            yin[2] = th2[i];
            yin[3] = w2[i];
         */

        double dydx[] = new double[N];
        /* function to fill array of derivatives dydx at xin */

        dydx[0] = yin[1];
        dydx[1] = -g/L1*Math.sin(th1);

        return dydx;

    }

    double[] runge_kutta(double xin, double yin[], double h)
    {
        double yout[] = new double[N];
        /* fourth order Runge-Kutta - see e.g. Numerical Recipes */

        int i;
        double hh, xh;
        double dydx[] = new double[N];
        double dydxt[] = new double[N];
        double yt[] = new double[N];
        double k1[] = new double[N];
        double k2[] = new double[N];
        double k3[] = new double[N];
        double k4[] = new double[N];

        hh = 0.5*h;
        xh = xin + hh;

        dydx = derivs(xin, yin); /* first step */
//        System.out.print("a:");
        for (i = 0; i < N; i++)
        {
            k1[i] = h*dydx[i];
            yt[i] = yin[i] + 0.5*k1[i];
//            System.out.printf(" %5.2f ",k1[i]);
        }

        dydxt = derivs(xh, yt); /* second step */
//        System.out.print("\nb:");
        for (i = 0; i < N; i++)
        {
            k2[i] = h*dydxt[i];
            yt[i] = yin[i] + 0.5*k2[i];
//            System.out.printf(" %5.2f ",k2[i]);
        } 

        dydxt = derivs(xh, yt); /* third step */
//        System.out.print("\nc:");
        for (i = 0; i < N; i++)
        {
            k3[i] = h*dydxt[i];
            yt[i] = yin[i] + k3[i];
//            System.out.printf(" %5.2f ",k3[i]);
        }

        dydxt = derivs(xin + h, yt); /* fourth step */
//        System.out.print("\nd:");
        for (i = 0; i < N; i++)
        {
            k4[i] = h*dydxt[i];
            yout[i] = yin[i] + k1[i]/6. + k2[i]/3. + k3[i]/3. + k4[i]/6.;
//            System.out.printf(" %5.2f ",k4[i]);
        }
//        System.out.println();
        return yout;

    }
}

class Pannello extends JPanel implements Runnable, MouseMotionListener{

    double x[];
    double y[];
    int offX, offY;
    double scale;
    final int DIM_BARRA_EN = 400;

    Sistema sistema;

    public Pannello(Sistema sistema){
        this.sistema=sistema;
        x = new double[2];
        y = new double[2];
    }

    @Override
    protected void paintComponent(Graphics g) {
        super.paintComponent(g);
        int raggioPallina=20;
        offX=this.getWidth()/2;
        offY=this.getHeight()/2;
        scale=(getHeight()-100)/sistema.L1/2;
        ArrayList<Point2D.Double> tracciato;
        synchronized(sistema){

            x[0]=sistema.getX()*scale;
            y[0]=-sistema.getY()*scale;
            tracciato = sistema.getTracciato();
        }        Graphics2D g2 = (Graphics2D)g;      
        g2.setStroke(new BasicStroke(3.0f));

        for(int i=0; i<tracciato.size()-1; i++){
            int x1 = (int) (tracciato.get(i).x*scale+offX);
            int y1 = (int) (-tracciato.get(i).y*scale+offY);
            int x2 = (int) (tracciato.get(i+1).x*scale+offX);
            int y2 = (int) (-tracciato.get(i+1).y*scale+offY);
            int c = (int) (255.0/tracciato.size()*i);
            //            System.out.println("i: "+i+" colore: "+c+" n."+sistema.tracciato.size()*i);
            g2.setPaint(new Color(c, 0, 0));
            g.drawLine(x1,y1,x2,y2);
        }

        g2.setPaint(new Color(0, 0, 0));
        g2.setStroke(new BasicStroke(2.0f));

        //bracci e palline
        g.drawLine(offX,offY, (int)x[0]+offX,(int)y[0]+offY);
        g.fillOval((int)x[0]-raggioPallina/2+offX, (int)y[0]-raggioPallina/2+offY, raggioPallina, raggioPallina);

        //energia
        g2.setStroke(new BasicStroke(1.0f));
        g2.setPaint(new Color(0.5f, 0.5f, 0.5f));
        g.fillRect(10,10,DIM_BARRA_EN,10);//tot
        g2.setPaint(new Color(0f, 1f, 0f));

        double e = sistema.energiaCinetica();
        double v = sistema.energiaPotenziale();
        double m = sistema.energiaIniziale;
        g.fillRect(10,20,(int)(e/sistema.energiaIniziale*DIM_BARRA_EN),10);//cinetica
        g2.setPaint(new Color(1f, 0f, 0f));
        g.fillRect(10+(int)(e/sistema.energiaIniziale*DIM_BARRA_EN),20,
                (int)(v/sistema.energiaIniziale*DIM_BARRA_EN),10);//potenziale    }

    @Override
    public void mouseDragged(MouseEvent arg0) {
        // TODO Auto-generated method stub

    }
    @Override
    public void mouseMoved(MouseEvent arg0) {

        x[0]=(double)arg0.getX()-offX;
        y[0]=(double)arg0.getY()-offY;

    }

    @Override
    public void run() {
        // TODO Auto-generated method stub
        while(true){
            sistema.simulate();

            System.out.format("EnIni:%8.2f Tot.%8.2f = C: %8.2f P: %8.2f\n",
                    sistema.energiaIniziale,
                    sistema.energiaCinetica()+sistema.energiaPotenziale(),
                    sistema.energiaCinetica(),sistema.energiaPotenziale());

            repaint();
            try {
                Thread.sleep((long) (sistema.dt*1000));
            } catch (InterruptedException e) {
                e.printStackTrace();
            }
        }
    }

}public class Pendolo {    public static void main(String[] args) {
        JFrame frame = new JFrame();
        frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        frame.setSize(600, 600);
        frame.setLayout(null);

        Sistema s = new Sistema();

        Pannello p = new Pannello(s);
        p.setBounds(0,0,600,600);
        frame.add(p);
        p.addMouseMotionListener(p);
        new Thread(p).start();

        s.simulate();

        frame.setVisible(true);
    }}

 

[lofgran]

Hi there,

Thank you for your question. It seems like you are experiencing a discrepancy between the results of your Runge-Kutta method and the results from the online pendulum simulation. There could be a few factors contributing to this difference.

Firstly, it is important to note that the online simulation may be using a different algorithm or method to simulate the pendulum. This could result in different results even with the same input values.

Additionally, it is possible that there may be a mistake in your implementation of the Runge-Kutta method. It may be helpful to compare your code with other sources or to use a software tool to check your calculations.

One potential software tool that you could use to compare your calculated values is MATLAB. It has built-in functions for solving differential equations, including the Runge-Kutta method. You could input your values into the MATLAB function and compare the results to your own.

I hope this helps and good luck with your studies!