みその計算物理学
ホーム はじめに リンク集
ルンゲクッタ法で様々な連立微分方程式を解く数値計算例(C言語)
  • プログラムソース (ラングフォードの方程式)
    /*ラングフォードの方程式
     dx/dt=(z-0.7)x-3.5y
      dy/dt=3.5x+(z-0.7)y
      dz/dt=a+z-z^3/3-(x^2+y^2)(1+0.25z)
      をルンゲクッタ法で解く
    */
    #include <stdio.h>
    #include <math.h>
    
    double f1(double t,double x,double y,double z);
    double f2(double t,double x,double y,double z);
    double f3(double t,double x,double y,double z);
    
    int main(void)
    {
            double t,x,y,z,dt,tmax;
            double k0[3],k1[3],k2[3],k3[3];
    
            dt=0.01;
            tmax=100.0;
            x=y=z=1.0;              //初期値
    
            FILE *output;
            output=fopen("output.data","w");
    
            for(t=0.0;t<tmax;t+=dt) {
                    k0[0]=dt*f1(t,x,y,z);
                    k0[1]=dt*f2(t,x,y,z);
                    k0[2]=dt*f3(t,x,y,z);
    
                    k1[0]=dt*f1(t+dt/2.0,x+k0[0]/2.0,y+k0[1]/2.0,z+k0[2]/2.0);
                    k1[1]=dt*f2(t+dt/2.0,x+k0[0]/2.0,y+k0[1]/2.0,z+k0[2]/2.0);
                    k1[2]=dt*f3(t+dt/2.0,x+k0[0]/2.0,y+k0[1]/2.0,z+k0[2]/2.0);
    
                    k2[0]=dt*f1(t+dt/2.0,x+k1[0]/2.0,y+k1[1]/2.0,z+k1[2]/2.0);
                    k2[1]=dt*f2(t+dt/2.0,x+k1[0]/2.0,y+k1[1]/2.0,z+k1[2]/2.0);
                    k2[2]=dt*f3(t+dt/2.0,x+k1[0]/2.0,y+k1[1]/2.0,z+k1[2]/2.0);
    
                    k3[0]=dt*f1(t+dt,x+k2[0],y+k2[1],z+k2[2]);
                    k3[1]=dt*f2(t+dt,x+k2[0],y+k2[1],z+k2[2]);
                    k3[2]=dt*f3(t+dt,x+k2[0],y+k2[1],z+k2[2]);
    
                    x=x+(k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0;
                    y=y+(k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0;
                    z=z+(k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0;
                    fprintf(output,"%f %f %f\n",x,y,z);
    
            }
    
            fclose(output);
    
            return 0;
    }
    
    double f1(double t,double x,double y,double z)
    {
            double r;
    
            r=(z-0.7)*x-3.5*y;
    
            return r;
    }
    
    double f2(double t,double x,double y,double z)
    {
            double r;
    
            r=3.5*x+(z-0.7)*y;
    
            return r;
    }
    
    double f3(double t,double x,double y,double z)
    {
            double r,a;
    
            a=0.6;          //この数値を変えることによって多様な解が得られる
            r=a+z-z*z*z/3.0-(x*x+y*y)*(1.0+0.25*z);
    
            return r;
    }
    
    
  • GNUPLOT出力結果
  • プログラムソース (ローレンツの方程式)
    /*ローレンツの方程式
      dx/dt=10(y-x)
      dy/dt=28x-y-xz
      dz/dt=-(8/3)z+xy
      をルンゲクッタ法で解く
    */
    #include <stdio.h>
    #include <math.h>
    
    double f1(double t,double x,double y,double z);
    double f2(double t,double x,double y,double z);
    double f3(double t,double x,double y,double z);
    
    int main(void)
    {
            double t,x,y,z,dt,tmax;
            double k0[3],k1[3],k2[3],k3[3];
    
            dt=0.01;
            tmax=100.0;
            x=y=z=1.0;              //初期値
    
            FILE *output;
            output=fopen("output.data","w");
    
            for(t=0.0;t<tmax;t+=dt) {
                    k0[0]=dt*f1(t,x,y,z);
                    k0[1]=dt*f2(t,x,y,z);
                    k0[2]=dt*f3(t,x,y,z);
    
                    k1[0]=dt*f1(t+dt/2.0,x+k0[0]/2.0,y+k0[1]/2.0,z+k0[2]/2.0);
                    k1[1]=dt*f2(t+dt/2.0,x+k0[0]/2.0,y+k0[1]/2.0,z+k0[2]/2.0);
                    k1[2]=dt*f3(t+dt/2.0,x+k0[0]/2.0,y+k0[1]/2.0,z+k0[2]/2.0);
    
                    k2[0]=dt*f1(t+dt/2.0,x+k1[0]/2.0,y+k1[1]/2.0,z+k1[2]/2.0);
                    k2[1]=dt*f2(t+dt/2.0,x+k1[0]/2.0,y+k1[1]/2.0,z+k1[2]/2.0);
                    k2[2]=dt*f3(t+dt/2.0,x+k1[0]/2.0,y+k1[1]/2.0,z+k1[2]/2.0);
    
                    k3[0]=dt*f1(t+dt,x+k2[0],y+k2[1],z+k2[2]);
                    k3[1]=dt*f2(t+dt,x+k2[0],y+k2[1],z+k2[2]);
                    k3[2]=dt*f3(t+dt,x+k2[0],y+k2[1],z+k2[2]);
    
                    x=x+(k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0;
                    y=y+(k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0;
                    z=z+(k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0;
                    fprintf(output,"%f %f %f\n",x,y,z);
    
            }
    
            fclose(output);
    
            return 0;
    }
    
    double f1(double t,double x,double y,double z)
    {
            double r;
    
            r=10*(y-x);
    
            return r;
    }
    
    double f2(double t,double x,double y,double z)
    {
            double r;
    
            r=28*x-y-x*z;
    
            return r;
    }
    
    double f3(double t,double x,double y,double z)
    {
            double r;
    
            r=-8/3*z+x*y;
    
            return r;
    }
    
    
  • GNUPLOT出力結果
  • プログラムソース (レスラーの方程式)
    /*レスラーの方程式
      dx/dt=-y-z
      dy/dt=x+0.34y
      dz/dt=0.4x-4.5z+xz
      をルンゲクッタ法で解く
    */
    #include <stdio.h>
    #include <math.h>
    
    double f1(double t,double x,double y,double z);
    double f2(double t,double x,double y,double z);
    double f3(double t,double x,double y,double z);
    
    int main(void)
    {
            double t,x,y,z,dt,tmax;
            double k0[3],k1[3],k2[3],k3[3];
    
            dt=0.01;
            tmax=300.0;
            x=y=z=1.0;              //初期値
    
            FILE *output;
            output=fopen("output.data","w");
    
            for(t=0.0;t<tmax;t+=dt) {
                    k0[0]=dt*f1(t,x,y,z);
                    k0[1]=dt*f2(t,x,y,z);
                    k0[2]=dt*f3(t,x,y,z);
    
                    k1[0]=dt*f1(t+dt/2.0,x+k0[0]/2.0,y+k0[1]/2.0,z+k0[2]/2.0);
                    k1[1]=dt*f2(t+dt/2.0,x+k0[0]/2.0,y+k0[1]/2.0,z+k0[2]/2.0);
                    k1[2]=dt*f3(t+dt/2.0,x+k0[0]/2.0,y+k0[1]/2.0,z+k0[2]/2.0);
    
                    k2[0]=dt*f1(t+dt/2.0,x+k1[0]/2.0,y+k1[1]/2.0,z+k1[2]/2.0);
                    k2[1]=dt*f2(t+dt/2.0,x+k1[0]/2.0,y+k1[1]/2.0,z+k1[2]/2.0);
                    k2[2]=dt*f3(t+dt/2.0,x+k1[0]/2.0,y+k1[1]/2.0,z+k1[2]/2.0);
    
                    k3[0]=dt*f1(t+dt,x+k2[0],y+k2[1],z+k2[2]);
                    k3[1]=dt*f2(t+dt,x+k2[0],y+k2[1],z+k2[2]);
                    k3[2]=dt*f3(t+dt,x+k2[0],y+k2[1],z+k2[2]);
    
                    x=x+(k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0;
                    y=y+(k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0;
                    z=z+(k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0;
                    fprintf(output,"%f %f %f\n",x,y,z);
    
            }
    
            fclose(output);
    
            return 0;
    }
    
    double f1(double t,double x,double y,double z)
    {
            double r;
    
            r=-y-z;
    
            return r;
    }
    
    double f2(double t,double x,double y,double z)
    {
            double r;
    
            r=x+0.34*y;
    
            return r;
    }
    
    double f3(double t,double x,double y,double z)
    {
            double r;
    
            r=0.4*x-4.5*z+x*z;
    
            return r;
    }
    
    
  • GNUPLOT出力結果