noname.c

偏微分方程四方法划曲线比较C源码. 欧拉法 改进欧拉法 龙格-库塔法 阿达姆斯法

#include <stdio.h>
#include <math.h>
#include <graphics.h>
#define a 0                /*****      区间     ******/
#define b 2                /**** 1.0—————6.4 ****/
#define h 0.1              /*区间步长0.002————0.1*/
#define num (int)((b-a)/h) /***** x区间分段总数 ******/
#define Yshang 0           /***** y坐标轴顶点坐标 ****/
#define Yxia 480           /***** y坐标轴下限坐标 ****/
#define Xzuo 0             /***** x坐标轴左坐标 ******/
#define Xyou 640           /***** x坐标轴右坐标 ******/
#define Euler_Color       14  /*****欧拉曲线颜色*****/
#define Adv_Euler_Color   11  /*****改进欧拉曲线颜色*****/
#define Runge_kutta_Color 10  /*****龙格库塔曲线颜色*****/
#define Adams_Color       12  /*****阿达姆斯曲线颜色*****/
#define mat_Color         15  /*****注释文字颜色*****/
#define mat_x 540             /*****注释文字坐标*****/
#define mat_y 30
#define Delaytime 5000        /*****画线延时大小*****/
/*#define PressAnyKey           /*****是否单步画线*****/*/
#define Delay                 /*****是否设置画线延时*****/

FILE *file;                /*****定义文件指针*****/
int Ymax,Xstep,Ystep,x;    /**** y最大约值、坐标轴x,y刻度长、图中x坐标值****/
int Euler_already=0;       /***** 判断欧拉法函数是否已运行的标志变量********/
int pixel_back=0;          /***** 是否将画笔撤回的标志变量**********/
/***************输出点阵字模*******************/
void drawmat(char *mat,int matsize,int x,int y,int color)
{int i,j,k,n;
 n=(matsize-1)/8+1;
 for(j=0;j<matsize;j++)
  for(i=0;i<n;i++)
   for(k=0;k<8;k++)
    if(mat[j*n+i]&(0x80>>k))
     putpixel(x+i*8+k,y+j,color);
}
/***************图表面板*******************/
void coordinate()
{ int i;
  /**********字模************/
  char mat_ou[24]= {0xF9,0x00,0x81,0x00,0x8B,0xE0,0xCA,0x20,
                    0xB5,0x00,0x91,0x00,0xA9,0x00,0xCA,0x80,
                    0x82,0x80,0xFC,0x40,0x08,0x20,0x00,0x00},
       mat_la[24]= {0x22,0x00,0x21,0x00,0xF7,0xE0,0x20,0x00,
                    0x24,0x40,0x34,0x40,0xE2,0x40,0x22,0x80,
                    0x22,0x80,0x21,0x00,0xEF,0xE0,0x00,0x00},
       mat_fa[24]= {0x41,0x00,0x21,0x00,0x2F,0xE0,0x81,0x00,
                    0x41,0x00,0x1F,0xE0,0x42,0x00,0x42,0x80,
                    0xC4,0x40,0x49,0xE0,0x4E,0x20,0x00,0x00},
       mat_gai[24]={0xF2,0x00,0x13,0xE0,0x12,0x40,0x14,0x40,
                    0xF4,0x40,0x8A,0x40,0x82,0x80,0x81,0x00,
                    0x99,0x80,0xE2,0x40,0x0C,0x20,0x00,0x00},
       mat_jin[24]={0x44,0x80,0x24,0x80,0x1F,0xE0,0x04,0x80,
                    0xC4,0x80,0x5F,0xE0,0x44,0x80,0x44,0x80,
                    0x48,0x80,0xB0,0x80,0x9F,0xE0,0x00,0x00},
       mat_a[24]  ={0xF7,0xE0,0x90,0x40,0xA0,0x40,0xC7,0x40,
                    0xA5,0x40,0x95,0x40,0x97,0x40,0xF0,0x40,
                    0x80,0x40,0x80,0x40,0x81,0xC0,0x00,0x00},
       mat_da[24] ={0x42,0x00,0x22,0x00,0x22,0x00,0x1F,0xE0,
                    0xC2,0x00,0x42,0x00,0x45,0x00,0x48,0xC0,
                    0x50,0x40,0xA0,0x00,0x9F,0xE0,0x00,0x00},
       mat_mu[24] ={0x27,0xC0,0x24,0x40,0xF5,0x40,0x55,0x40,
                    0x5F,0xE0,0x54,0x40,0x55,0x40,0xD5,0x40,
                    0x27,0xE0,0x50,0x40,0x81,0x80,0x00,0x00},
       mat_si[24] ={0x48,0x60,0xFD,0x80,0x49,0x00,0x79,0x00,
                    0x49,0xE0,0x79,0x40,0x49,0x40,0xFD,0x40,
                    0x2A,0x40,0x4A,0x40,0x84,0x40,0x00,0x00},
       mat_long[24]={0x09,0x00,0x08,0x80,0xFF,0xE0,0x0A,0x00,
                    0x0A,0x40,0x0A,0x80,0x13,0x00,0x12,0x00,
                    0x26,0x20,0x5A,0x20,0x83,0xE0,0x00,0x00},
       mat_ge[24] ={0x22,0x00,0x23,0xC0,0xF4,0x40,0x2A,0x80,
                    0x61,0x00,0x72,0x80,0xAC,0x60,0xA7,0xC0,
                    0x24,0x40,0x24,0x40,0x27,0xC0,0x00,0x00},
       mat_ku[24] ={0x02,0x00,0x7F,0xE0,0x44,0x00,0x7F,0xE0,
                    0x4A,0x00,0x4A,0x00,0x5F,0xC0,0x42,0x00,
                    0x7F,0xE0,0x82,0x00,0x82,0x00,0x00,0x00},
       mat_ta[24] ={0x44,0x80,0x5F,0xE0,0x44,0x80,0xF3,0x00,
                    0x44,0x80,0x4F,0xC0,0x50,0x20,0x47,0xC0,
                    0x74,0x40,0xC4,0x40,0x07,0xC0,0x00,0x00};

  setbkcolor(1);                      /********设置图表背景色*********/
  line(Xzuo+40,Yshang,Xzuo+40,Yxia);  /*****绘制坐标轴及其方向********/
  line(Xzuo,Yxia-40,Xyou,Yxia-40);
  moveto(Xzuo+40,Yshang);
  lineto(Xzuo+43,Yshang+15);
  lineto(Xzuo+40,Yshang+10);
  lineto(Xzuo+37,Yshang+15);
  lineto(Xzuo+40,Yshang);
  moveto(Xyou,Yxia-40);
  lineto(Xyou-17,Yxia-43);
  lineto(Xyou-10,Yxia-40);
  lineto(Xyou-17,Yxia-37);
  lineto(Xyou,Yxia-40);
  settextstyle(3,0,1);
  outtextxy(Xyou-30,Yxia-30,"X");
  outtextxy(Xzuo+20,Yshang+25,"Y");
  outtextxy(Xzuo+20,Yxia-30,"O");
  Xstep=((Xyou-40)-(Xzuo+40))/num;            /********计算x,y轴刻度长**********/
  Ystep=((Yxia-40)-(Yshang+200))/Ymax;
  for(i=Xzuo+40+Xstep;i<Xyou-40;i+=Xstep)    /********绘制x轴刻度*********/
  { line(i,Yxia-40,i,Yxia-38);
    }
  for(i=Yxia-40-Ystep;i>=Yshang+200;i-=Ystep) /********绘制y轴刻度*********/
  { line(Xzuo+40,i,Xzuo+38,i);
    }
  drawmat(mat_ou,12,mat_x,mat_y,mat_Color);       /****输出字模“欧拉法”*********/
  drawmat(mat_la,12,mat_x+13,mat_y,mat_Color);
  drawmat(mat_fa,12,mat_x+26,mat_y,mat_Color);

  drawmat(mat_gai,12,mat_x,mat_y+15,mat_Color);   /****输出字模“改进欧拉法”*****/
  drawmat(mat_jin,12,mat_x+13,mat_y+15,mat_Color);
  drawmat(mat_ou,12,mat_x+26,mat_y+15,mat_Color);
  drawmat(mat_la,12,mat_x+39,mat_y+15,mat_Color);
  drawmat(mat_fa,12,mat_x+52,mat_y+15,mat_Color);

  drawmat(mat_long,12,mat_x,mat_y+30,mat_Color);  /****输出字模“龙格库塔法”*****/
  drawmat(mat_ge,12,mat_x+13,mat_y+30,mat_Color);
  drawmat(mat_ku,12,mat_x+26,mat_y+30,mat_Color);
  drawmat(mat_ta,12,mat_x+39,mat_y+30,mat_Color);
  drawmat(mat_fa,12,mat_x+52,mat_y+30,mat_Color);

  drawmat(mat_a,12,mat_x,mat_y+45,mat_Color);     /****输出字模“阿达姆斯法”*****/
  drawmat(mat_da,12,mat_x+13,mat_y+45,mat_Color);
  drawmat(mat_mu,12,mat_x+26,mat_y+45,mat_Color);
  drawmat(mat_si,12,mat_x+39,mat_y+45,mat_Color);
  drawmat(mat_fa,12,mat_x+52,mat_y+45,mat_Color);

  setcolor(Euler_Color);       line(mat_x-7,mat_y+6,mat_x-70,mat_y+6);  /**绘制“欧拉法”标示线**/
  setcolor(Adv_Euler_Color);   line(mat_x-7,mat_y+21,mat_x-70,mat_y+21);/**绘制“改进欧拉法”标示线**/
  setcolor(Runge_kutta_Color); line(mat_x-7,mat_y+36,mat_x-70,mat_y+36);/**绘制“龙格库塔法”标示线**/
  setcolor(Adams_Color);       line(mat_x-7,mat_y+51,mat_x-70,mat_y+51);/**绘制“阿达姆斯法”标示线**/
 }
/***************绘制曲线*******************/
void draw(double u,int color)
{ int y;
  y=(Yxia-40)-(int)(u*Ystep);       /**********将函数值u转换成坐标值y************/
  if(pixel_back==1)    /*******画笔撤回，将x坐标值重置为x轴原点********/
  { x=Xzuo+40;
    putpixel(x,y,color);
    setcolor(color);
    moveto(x,y);
    pixel_back=0;
    }
  else
    lineto(x,y);
  x=x+Xstep;
 }
/***************函数表达式***************/
double f(double t,double u)
{return(exp(t)-2*u);
 }
/***************欧拉法*******************/
double Euler(double t0,double u0)
{double i,u1=u0,t1=t0,u2,t2;
 int n=1,color=Euler_Color;
 pixel_back=1;
 if(Euler_already==1)
 { fprintf(file,"/***************欧拉法*******************/\n");
  }
 for(i=a;i<=b;i+=h)
 { u2=u1+h*f(t1,u1);          /**********“欧拉法”计算公式**********/
   t2=t1+h;
/********首次调用(Euler_already==0)此函数时只返回y的最大约值，
   再次调用(Euler_already==1)时输出计算结果并绘图**********/
   if(Euler_already==1)
   { fprintf(file,"t%d=%lf, u%d=%lf\n",n,t2,n,u2);
     #ifdef Delay
      delay(Delaytime);       /*********延时绘图条件编译***********/
     #endif
     draw(u2,color);          /**********立即画线***********/
    }
   n++;
   u1=u2;
   t1=t2;
   }
   Euler_already=1;
   return(u1);

}
/***************改进欧拉法*******************/
void Adv_Euler(double t0,double u0)
{double i,u1=u0,t1=t0,u2,t2,u3,t3=t0;
 int n=1,color=Adv_Euler_Color;
 pixel_back=1;
 fprintf(file,"/***************改进欧拉法*******************/\n");
 u2=u1+h*f(t1,u1);                   /********利用“欧拉法”计算初值************/
 t2=t1+h;
 for(i=a;i<=b;i+=h)
 { u3=u2+(h/2)*(f(t1,u1)+f(t2,u2));  /***********“改进欧拉法”计算公式*********/
   t3=t3+h;
   fprintf(file,"t%d=%lf, u%d=%lf\n",n,t3,n,u3);
     #ifdef Delay
      delay(Delaytime);              /*********延时绘图条件编译***********/
     #endif
   draw(u3,color);                   /**********立即画线***********/
   n++;
   u1=u2;t1=t2;
   u2=u3;t2=t3;
   }
}
/***************龙格-库塔法*******************/
void Runge_Kutta(double t0,double u0)
{int n=1,color=Runge_kutta_Color;
 double i;
 double k1,k2,k3,k4;
 double t1=t0,u1=u0,t2,u2;
 pixel_back=1;
 fprintf(file,"/***************龙格-库塔法*******************/\n");
 for(i=a;i<=b;i+=h)
 { k1=h*f(t1,u1);         /*************“龙格-库塔法”计算公式*************/
   k2=h*f(t1+h/2,u1+k1/2);
   k3=h*f(t1+h/2,u1+k2/2);
   k4=h*f(t1+h,u1+k3);
   u2=u1+(k1+(2*k2)+(2*k3)+k4)/6.0;
   t2=t1+h;
   fprintf(file,"t%d=%lf, u%d=%lf\n",n,t2,n,u2);
     #ifdef Delay
      delay(Delaytime);   /*********延时绘图条件编译***********/
     #endif
   draw(u2,color);        /**********立即画线***********/
   n++;
   u1=u2;t1=t2;
   }
}
/***************阿达姆斯法*******************/
void Adams(double t0,double u0)
{ double i,u1=u0,t1=t0,u2,t2,u3,t3=t0;
  double k1,k2,k3,k4;
  int n=1,color=Adams_Color;
  pixel_back=1;
  k1=h*f(t1,u1);         /*********利用“龙格库塔法”计算初值*********/
  k2=h*f(t1+h/2,u1+k1/2);
  k3=h*f(t1+h/2,u1+k2/2);
  k4=h*f(t1+h,u1+k3);
  u2=u1+(k1+(2*k2)+(2*k3)+k4)/6.0;
  t2=t1+h;
  fprintf(file,"/***************阿达姆斯法*******************/\n");
 for(i=a;i<=b;i+=h)
 { u3=u2+(h/2)*(3*f(t2,u2)-f(t1,u1)); /**********“阿达姆斯法”计算公式*************/
   t3=t3+h;
   fprintf(file,"t%d=%lf, u%d=%lf\n",n,t3,n,u3);
     #ifdef Delay
      delay(Delaytime);               /*********延时绘图条件编译***********/
     #endif
   draw(u3,color);                    /**********立即画线***********/
   n++;
   u1=u2;t1=t2;
   u2=u3;t2=t3;
   }
}
/***************主函数******************/
main()
{double x0,y0;
 int gd=0,gm=1;
 initgraph(&gd,&gm,"");     /*********初始化图形模式***********/
 x0=0;
 y0=1;
 file=fopen("result.txt","w"); /**将运算结果输出到result.txt文本中**/
 Ymax=Euler(x0,y0)+1;          /*******取得y最大约值******/
 coordinate();              /**********绘制图表板面************/
 #ifdef PressAnyKey
   getch();                    /*****单步画线条件编译*****/
 #endif
 Euler(x0,y0);                 /*********“欧拉法”***********/
 #ifdef PressAnyKey
   getch();
 #endif
 Adv_Euler(x0,y0);             /*********“改进欧拉法”*******/
 #ifdef PressAnyKey
   getch();
 #endif
 Runge_Kutta(x0,y0);           /*********“龙格库塔法”*******/
 #ifdef PressAnyKey
   getch();
 #endif
 Adams(x0,y0);                 /*********“阿达姆斯法”*******/
 fclose(file);                 /***********关闭文件***************/
 getch();
 closegraph();                 /*********关闭图形模式***************/
}