runge-kuttac++.txt

在使用龙格-库塔(RK)方法对连续系统进行数字仿真时,为了保证数值计算的稳定性以及仿真结果具有足够的精度,通常采用变步长策略。为了有效地解决变步长仿真计算过程中,输出节点与计算节点不相吻合的问题,该文

隐式Runge-Kutta的C++程序



#define N 21
#include "imrk63.c"
void f( double x, double y, double ffd[2] );
extern void imrk63( int iter, double x0, double x1, double y[N] );
void main( )
{
int i,iter;
double x0,x1,h,x,y0,y[N];
iter=5;
x0=0.0;  x1=1.0;
y0=0.0;  y[0]=y0;
system("cls");
printf("Sixth-order implicit Runge-Kutta method:\n");
imrk63( iter, x0, x1, y);
h=(x1-x0)/(N-1.0);
printf("     xn           yn\n");
for(i=0;i<N;i++)
  { x=x0+i*h; printf("%10.6f %12.8f\n",x,y); }
getch();
}
/********* function f( x, y, ffd ) ********/
void f( double x, double y, double ffd[2] )
{
ffd[0]=y*y+x*x;
ffd[1]=2.0*y;
} 
 
2007-9-7 17:49 #1 
     
  
chunming2005 

工程师 




精华 0
积分 57
帖子 95
水位 190 
技术分 0 
状态 离线  复制内容到剪贴板代码:/*imrk63.c: Sixth-order third stage implicit Runge-kutta method. */
#define DIM 3
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int gcpelim( double A[DIM][DIM], double xx[DIM] );
extern void f( double x, double y, double ffd[2] );
void imrk63( int iter, double x0, double x1, double y[N] )
{
int i,j,ij,ik;
double xx,yy,k[3],h,ac[3],bc[3][3],c[3],A[3][3];
double ff[3],fd[3],ffd[2];
ac[0]=(5.0-sqrt(15.0))/10.0;
ac[1]=0.5;
ac[2]=(5.0+sqrt(15.0))/10.0;
bc[0][0]=5.0/36.0;
bc[0][1]=(10.0-3.0*sqrt(15.0))/45.0;
bc[0][2]=(25.0-6.0*sqrt(15.0))/180.0;
bc[1][0]=(10.0+3.0*sqrt(15.0))/72.0;
bc[1][1]=2.0/9.0;
bc[1][2]=(10.0-3.0*sqrt(15.0))/72.0;
bc[2][0]=(25.0+6.0*sqrt(15.0))/180.0;
bc[2][1]=(10.0+3.0*sqrt(15.0))/45.0;
bc[2][2]=5.0/36.0;
k[0]=0.0;  k[1]=0.0;  k[2]=0.0;
h=(x1-x0)/(N-1.0);
for(i=0;i<N-1;i++)
  {
  for(ik=0;ik<iter;ik++)
    {
    for(j=0;j<3;j++)
      {
      xx=x0+i*h+ac[j]*h;
      yy=y+bc[j][0]*k[0]+bc[j][1]*k[1]+bc[j][2]*k[2];
      f( xx, yy, ffd );
      ff[j]=h*ffd[0];
      fd[j]=h*ffd[1];
      }
    for(ij=0;ij<3;ij++)
      {
      for(j=0;j<3;j++)  A[ij][j]=-fd[ij]*bc[ij][j];
      A[ij][ij]=1.0+A[ij][ij];
      }
    for(ij=0;ij<3;ij++)  c[ij]=ff[ij]-k[ij];
    if( gcpelim( A, c ) == 1 )
      {
      printf("The linear system hasn't solution!\n");
      printf("Strike any key to exit!\n"); getch(); exit(1);
      }
    for(j=0;j<3;j++) k[j]=c[j]+k[j];
    }
  y[i+1]=y+(5.0*k[0]+8.0*k[1]+5.0*k[2])/18.0;
  }
}
/****** Gauss column principal eliminate for linear system ********/
#define EPSILON 0.000000001
int gcpelim( double A[DIM][DIM], double xx[DIM] )
{
int k,i,j,i0;
double pelement;
for(k=0;k<DIM;k++)
  {
  pelement=fabs(A[k][k]);       i0=k;
  for(i=k;i<DIM;i++)
    if( fabs(A[k]) > pelement ) { pelement=fabs(A[k]); i0=i; }
  if( i0 != k )
    {
    for(j=0;j<DIM;j++)
      { pelement=A[k][j]; A[k][j]=A[i0][j]; A[i0][j]=pelement; }
    pelement=xx[k]; xx[k]=xx[i0]; xx[i0]=pelement;
    }
  if( fabs(A[k][k]) < EPSILON ) return(1);
  for(i=k+1;i<DIM;i++)
    {
    A[k]=A[k]/A[k][k];
    for(j=k+1;j<DIM;j++) A[j]=A[j]-A[k]*A[k][j];
    xx=xx-A[k]*xx[k];
    }
  }
for(i=DIM-1;i>=0;i--)
  {
  for(j=i+1;j<DIM;j++) xx=xx-A[j]*xx[j];
  xx=xx/A;
  }
return(0);
}