Runge-Kutta method for Simple pendulum in java

Question

我正在尝试研究龙格-库塔方法并应用于简单的钟摆。使用时间步长 dt=0.1 (h=0.1(，钟摆每周期增加 1% 的能量......在本站点中：链接能量降低约0.01% 我做错了什么？有一个软件可以在哪里比较我的计算值？谢谢

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.2fn",
                    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();

        frame.setVisible(true);
    }

}

参考

Answer 1

要根据 RK4 理论获得结果，您应该使用参数中给出的状态来计算导数，而不是最后一个采样点的模型状态的全局变量。

double[] derivs(double xin, double yin[])
{
    double dydx[] = new double[N];
    /* function to fill array of derivatives dydx at xin */
    dydx[0] = yin[1]; 
    dydx[1] = -g/L1*Math.sin(yin[0]);
    /*        Error was here ^^^^^^^     */
    return dydx;
}

---

Runge-Kutta 方法适用于f(x,y) y'(x)=f(x,y(x))函数，这些函数只应完全依赖于这些输入变量或对整个积分全局的常量。

您之前的代码在每次derivs评估中引入了顺序h错误，因为您使用y[0]而不是y[0]+0.5*k1[0]等，这在 RK 步骤级别转换为二次误差。这意味着该方法与显式欧拉方法一样好，全局误差顺序为 1。