隐式Runge-Kutta的C++程序
/*imrk63m.c: Sixth-order third stage implicit Runge-Kutta method. */#define N 21
#include "imrk63.c"
void f( double x, double y, double ffd );
extern void imrk63( int iter, double x0, double x1, double y );
void main( )
{
int i,iter;
double x0,x1,h,x,y0,y;
iter=5;
x0=0.0;x1=1.0;
y0=0.0;y=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 )
{
ffd=y*y+x*x;
ffd=2.0*y;
} /*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, double xx );
extern void f( double x, double y, double ffd );
void imrk63( int iter, double x0, double x1, double y )
{
int i,j,ij,ik;
double xx,yy,k,h,ac,bc,c,A;
double ff,fd,ffd;
ac=(5.0-sqrt(15.0))/10.0;
ac=0.5;
ac=(5.0+sqrt(15.0))/10.0;
bc=5.0/36.0;
bc=(10.0-3.0*sqrt(15.0))/45.0;
bc=(25.0-6.0*sqrt(15.0))/180.0;
bc=(10.0+3.0*sqrt(15.0))/72.0;
bc=2.0/9.0;
bc=(10.0-3.0*sqrt(15.0))/72.0;
bc=(25.0+6.0*sqrt(15.0))/180.0;
bc=(10.0+3.0*sqrt(15.0))/45.0;
bc=5.0/36.0;
k=0.0;k=0.0;k=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*h;
yy=y+bc*k+bc*k+bc*k;
f( xx, yy, ffd );
ff=h*ffd;
fd=h*ffd;
}
for(ij=0;ij<3;ij++)
{
for(j=0;j<3;j++)A=-fd*bc;
A=1.0+A;
}
for(ij=0;ij<3;ij++)c=ff-k;
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=c+k;
}
y=y+(5.0*k+8.0*k+5.0*k)/18.0;
}
}
/****** Gauss column principal eliminate for linear system ********/
#define EPSILON 0.000000001
int gcpelim( double A, double xx )
{
int k,i,j,i0;
double pelement;
for(k=0;k<DIM;k++)
{
pelement=fabs(A); i0=k;
for(i=k;i<DIM;i++)
if( fabs(A) > pelement ) { pelement=fabs(A); i0=i; }
if( i0 != k )
{
for(j=0;j<DIM;j++)
{ pelement=A; A=A; A=pelement; }
pelement=xx; xx=xx; xx=pelement;
}
if( fabs(A) < EPSILON ) return(1);
for(i=k+1;i<DIM;i++)
{
A=A/A;
for(j=k+1;j<DIM;j++) A=A-A*A;
xx=xx-A*xx;
}
}
for(i=DIM-1;i>=0;i--)
{
for(j=i+1;j<DIM;j++) xx=xx-A*xx;
xx=xx/A;
}
return(0);
}
本程序来自于： 科学计算和C程序集 谢谢楼主分享
页:
[1]