隐式Runge-Kutta的C++程序

风花雪月

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. }

复制代码

风花雪月

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程序集

sizhiyuan2006

谢谢楼主分享