/*
continuous time lag equation

dx/dt = -a*x(t-r), t\geq 0
x(0)=1.0

r=8
r=10
*/

#include 
#include 
#include 

#define h 0.1
#define r 10 /*time lag*/
#define a 2.0 /*parameter*/

double FUNC();

main()
{
  int j; /*j; step number*/
  double nx[151],x[151];
  double nt=0,t=0.0;
  double k1=0.0,k2=0.0,k3=0.0,k4=0.0;

  for(j=0;j<150;j++)
    {
      if(j<=r){
	x[j]=1.0;
	/*printf("%d %f\n",j, x[j]);*/
      }
      else{
       
	k1=h*FUNC(t,x[j-r]);
	k2=h*FUNC(t+0.5*h,x[j-r]+0.5*k1);
	k3=h*FUNC(t+0.5*h,x[j-r]+0.5*k2);
	k4=h*FUNC(t+h,x[j-r]+k3);
	
	nx[j+1]=x[j]+0.166666667*(k1+2.0*k2+2.0*k3+k4);
	nt=t+h;
      /*printf("%f\n",x[j]);*/
      printf("%d %f\n",j-r, x[j]);
      x[j+1]=nx[j+1];
      t=nt;
      }
    }
  return 0;
}

/*FUNCTION */
double FUNC(double t, double x)
{
  return ( -a*x );
}


出力結果

RED : r=0(no time lag), GREEN : Oscillation , BLUE : Divergence


Back to C Language

Google




BLOG
PICASAWEB
Panoramio


REF.


タイムラグをもつ微分方程式 - 関数微分方程式入門