単純な振り子のRunge-Kuttaメソッド

Question

私はRunge-Kuttaメソッドを学び、単純な振り子に適用しようとしています。タイムステップdt = 0.1（h = 0.1）を使用すると、振り子は毎周期1％のエネルギーを発生します。このサイトでは： link エネルギーを約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.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(); 



     frame.setVisible(true); 
    } 


} 

Reference

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; 
} 

---

ルンゲクッタ法は、正確にこれらの入力変数にまたは全体の統合へのグローバルな定数に依存すべきy'(x)=f(x,y(x))で機能f(x,y)のために作られています。かわりにRK段階のレベルで二次誤差に変換等y[0]+0.5*k1[0]のy[0]を使用するように

以前コードはderivsの評価のそれぞれにおける順序hの誤差を導入します。これは、メソッドが明示的なEulerメソッドと同じくらい良いことを意味します。グローバルエラーオーダー1を持ちます。