查看: 2212|回复: 2

[C/C++] 隐式Runge-Kutta的C++程序

[复制链接]
发表于 2006-11-12 06:07 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
  1. /*imrk63m.c: Sixth-order third stage implicit Runge-Kutta method. */
  2. #define N 21
  3. #include "imrk63.c"
  4. void f( double x, double y, double ffd[2] );
  5. extern void imrk63( int iter, double x0, double x1, double y[N] );
  6. void main( )
  7. {
  8. int i,iter;
  9. double x0,x1,h,x,y0,y[N];
  10. iter=5;
  11. x0=0.0;  x1=1.0;
  12. y0=0.0;  y[0]=y0;
  13. system("cls");
  14. printf("Sixth-order implicit Runge-Kutta method:\n");
  15. imrk63( iter, x0, x1, y);
  16. h=(x1-x0)/(N-1.0);
  17. printf("     xn           yn\n");
  18. for(i=0;i<N;i++)
  19.   { x=x0+i*h; printf("%10.6f %12.8f\n",x,y[i]); }
  20. getch();
  21. }
  22. /********* function f( x, y, ffd ) ********/
  23. void f( double x, double y, double ffd[2] )
  24. {
  25. ffd[0]=y*y+x*x;
  26. ffd[1]=2.0*y;
  27. }
复制代码
回复
分享到:

使用道具 举报

 楼主| 发表于 2006-11-12 06:07 | 显示全部楼层
  1. /*imrk63.c: Sixth-order third stage implicit Runge-kutta method. */
  2. #define DIM 3
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #include <math.h>
  6. int gcpelim( double A[DIM][DIM], double xx[DIM] );
  7. extern void f( double x, double y, double ffd[2] );
  8. void imrk63( int iter, double x0, double x1, double y[N] )
  9. {
  10. int i,j,ij,ik;
  11. double xx,yy,k[3],h,ac[3],bc[3][3],c[3],A[3][3];
  12. double ff[3],fd[3],ffd[2];
  13. ac[0]=(5.0-sqrt(15.0))/10.0;
  14. ac[1]=0.5;
  15. ac[2]=(5.0+sqrt(15.0))/10.0;
  16. bc[0][0]=5.0/36.0;
  17. bc[0][1]=(10.0-3.0*sqrt(15.0))/45.0;
  18. bc[0][2]=(25.0-6.0*sqrt(15.0))/180.0;
  19. bc[1][0]=(10.0+3.0*sqrt(15.0))/72.0;
  20. bc[1][1]=2.0/9.0;
  21. bc[1][2]=(10.0-3.0*sqrt(15.0))/72.0;
  22. bc[2][0]=(25.0+6.0*sqrt(15.0))/180.0;
  23. bc[2][1]=(10.0+3.0*sqrt(15.0))/45.0;
  24. bc[2][2]=5.0/36.0;
  25. k[0]=0.0;  k[1]=0.0;  k[2]=0.0;
  26. h=(x1-x0)/(N-1.0);
  27. for(i=0;i<N-1;i++)
  28.   {
  29.   for(ik=0;ik<iter;ik++)
  30.     {
  31.     for(j=0;j<3;j++)
  32.       {
  33.       xx=x0+i*h+ac[j]*h;
  34.       yy=y[i]+bc[j][0]*k[0]+bc[j][1]*k[1]+bc[j][2]*k[2];
  35.       f( xx, yy, ffd );
  36.       ff[j]=h*ffd[0];
  37.       fd[j]=h*ffd[1];
  38.       }
  39.     for(ij=0;ij<3;ij++)
  40.       {
  41.       for(j=0;j<3;j++)  A[ij][j]=-fd[ij]*bc[ij][j];
  42.       A[ij][ij]=1.0+A[ij][ij];
  43.       }
  44.     for(ij=0;ij<3;ij++)  c[ij]=ff[ij]-k[ij];
  45.     if( gcpelim( A, c ) == 1 )
  46.       {
  47.       printf("The linear system hasn't solution!\n");
  48.       printf("Strike any key to exit!\n"); getch(); exit(1);
  49.       }
  50.     for(j=0;j<3;j++) k[j]=c[j]+k[j];
  51.     }
  52.   y[i+1]=y[i]+(5.0*k[0]+8.0*k[1]+5.0*k[2])/18.0;
  53.   }
  54. }
  55. /****** Gauss column principal eliminate for linear system ********/
  56. #define EPSILON 0.000000001
  57. int gcpelim( double A[DIM][DIM], double xx[DIM] )
  58. {
  59. int k,i,j,i0;
  60. double pelement;
  61. for(k=0;k<DIM;k++)
  62.   {
  63.   pelement=fabs(A[k][k]);       i0=k;
  64.   for(i=k;i<DIM;i++)
  65.     if( fabs(A[i][k]) > pelement ) { pelement=fabs(A[i][k]); i0=i; }
  66.   if( i0 != k )
  67.     {
  68.     for(j=0;j<DIM;j++)
  69.       { pelement=A[k][j]; A[k][j]=A[i0][j]; A[i0][j]=pelement; }
  70.     pelement=xx[k]; xx[k]=xx[i0]; xx[i0]=pelement;
  71.     }
  72.   if( fabs(A[k][k]) < EPSILON ) return(1);
  73.   for(i=k+1;i<DIM;i++)
  74.     {
  75.     A[i][k]=A[i][k]/A[k][k];
  76.     for(j=k+1;j<DIM;j++) A[i][j]=A[i][j]-A[i][k]*A[k][j];
  77.     xx[i]=xx[i]-A[i][k]*xx[k];
  78.     }
  79.   }
  80. for(i=DIM-1;i>=0;i--)
  81.   {
  82.   for(j=i+1;j<DIM;j++) xx[i]=xx[i]-A[i][j]*xx[j];
  83.   xx[i]=xx[i]/A[i][i];
  84.   }
  85. return(0);
  86. }
复制代码


本程序来自于: 科学计算和C程序集
发表于 2017-10-9 13:45 | 显示全部楼层
谢谢楼主分享
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表