みその計算物理学
ホーム はじめに リンク集
エムデン方程式の数値計算例(Java)
  • プログラム例
    import java.io.*;
    
    /**
     * エムデン方程式を解くプログラム。
     */
    
    public class Emdenequation
    {
            static double n;
    
            public static void main()
            {
            double t,t0,x,v,dt,tmax;
            double k1[]=new double[2],k2[]=new double[2],k3[]=new double[2],k4[]=new double[2];
                    String s;
    
                    Emdenequation f = new Emdenequation();
    
                    System.out.println("ポリトロープ指数の値は?");
                    try
                    {
                            InputStreamReader isr = new InputStreamReader(System.in);
                            BufferedReader br = new BufferedReader(isr);
                            s = br.readLine();
                            n = Double.parseDouble(s);
                            br.close();
                    }
                    catch (Exception e)
                    {
                            System.out.println(e);
                    }
    
                    t0 = 0.0001;      //初期値
                    x = 1.0 - 1.0 / 6.0 * t0 * t0 + n / 120.0 * t0 * t0 * t0 * t0;        //初期値
                    v = -1.0 / 3.0 * t0 + n / 30.0 * t0 * t0 * t0;  //初期値
                    dt = 0.001;       //刻み幅
                    tmax = 10000.0;   //最大値
    
                    try
                    {
                            FileWriter fw = new FileWriter("output.data");  //gnuplot表示用データファイル
    
                            for (t = t0; t < tmax; t += dt)
                            {
                                    k1[0] = dt * f.f1(t, x, v);
                                    k1[1] = dt * f.f2(t, x, v);
    
                                    k2[0] = dt * f.f1(t + dt / 2.0, x + k1[0] / 2.0, v + k1[1] / 2.0);
                                    k2[1] = dt * f.f2(t + dt / 2.0, x + k1[0] / 2.0, v + k1[1] / 2.0);
    
                                    k3[0] = dt * f.f1(t + dt / 2.0, x + k2[0] / 2.0, v + k2[1] / 2.0);
                                    k3[1] = dt * f.f2(t + dt / 2.0, x + k2[0] / 2.0, v + k2[1] / 2.0);
    
                                    k4[0] = dt * f.f1(t + dt, x + k3[0], v + k3[1]);
                                    k4[1] = dt * f.f2(t + dt, x + k3[0], v + k3[1]);
    
                                    x = x + (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0;
                                    v = v + (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0;
    
                                    fw.write(t + " " + x + " " + v + "\n");
    
    
                                    if (x < 0) break;
                            }
    
                    }
                    catch(Exception e)
                    {
                            System.out.print(e);
                    }
            }
            double f1(double t, double x, double v)
            {
                    return v;
            }
    
            double f2(double t, double x, double v)
            {
                    return -Math.pow(x, n) - 2.0 / t * v;
            }
    }
    
  • 数値計算結果

    以下に上記のプログラムで求めた解を示す。横軸はξで縦軸はθである。ポリトロピック指数は0とした。