****************************
* CMP203 Spring 98 *
* Homework # 1 *
* *
* Qin Shen *
* jeanshen@cse *
****************************
>>>> C version <<<<
#define INVALID_NUM -100000
#define EPSILON 0.000001
#define N 3
#define true 1
#define false 0
double x, xE;
double y[N], yE[N], z[N];
/* sample functions defined by me, the return value = f(x,y) = y' */
double FUNC( int which, double x, double y )
{
switch ( which )
{
case 0: return( 2*x );
case 1: return( 3*x*x );
case 2: return( y );
default:
printf("No this function.\n");
return( INVALID_NUM );
}
}
void FKT( double x, double *y, int n, double *z )
{
int i;
for ( i = 0; i < n; i++ )
z[i] = FUNC( i, x, y[i] );
}
void RK1ST( double x, double *y, double h, int n, double *xe, double *ye )
{
int j, k;
double w[N];
double a[5];
a[0] = a[1] = a[4] = h/2; a[2] = a[3] = h;
*xe = x;
for ( k = 0; k < n; k++ )
ye[k] = w[k] = y[k];
for ( j = 0; j < 5; j++ )
{
FKT( *xe, w, n, z );
*xe = x + a[j];
for( k = 0; k < n; k++ )
{
w[k] = y[k] + a[j] * z[k];
ye[k] = ye[k] + a[j+1] * z[k] / 3;
}
}
}
void RK( double x, double *y, int n, double eps, double xE,
double *yE, int fi )
{
double y1[N], y2[N], y3[N];
double *x1, *x2, *x3, H;
static double s, Hs;
int out, j, k;
int do_BB, do_EE; /* dummy flags - corresponding to Algol60 code */
if ( fi )
{
H = xE - x;
s = 0;
}
else
{
H = Hs;
out = false;
}
/* AA */
while ( !out )
{
if ( (( x + 2.01 * H - xE ) > 0) == ( H > 0) )
{
Hs = H;
out = true;
H = ( xE - x ) / 2;
}
RK1ST( x, y, 2*H, n, x1, y1 );
do_BB = true;
/* BB */
while ( do_BB )
{
do_EE = true;
RK1ST( x, y, H, n, x2, y2 );
RK1ST( *x2, y2, H, n, x3, y3 );
for ( k = 0; k < n; k++ )
{
if ( ( y1[k] - y3[k] ) > eps )
{
/* CC */
H = 0.5 * H;
out = false;
*x1 = *x2;
for ( j = 0; j < n; j++ )
y1[j] = y2[j];
do_EE = false;
break; /* goto BB */
}
}
/* EE - second part of BB */
if ( do_EE )
{
x = *x3;
if ( out )
{
for ( k = 0; k < n; k++ )
yE[k] = y3[k];
}
else
{
for ( k = 0; k < n; k++ )
y[k] = y3[k];
if ( s == 5 )
{
s = 0;
H = 2 * H;
}
s++;
do_BB = false;
}
}
}
}
}
int main( )
{
/* should initialize here */
RK( x, y, N, EPSILON, xE, yE, true );
/* should print information here */
}
--------------------------------------------------
>>>> Pascal version <<<<
program Runge_Kutta( input, output );
const INVALID_NUM = -100000;
EPSILON = 0.000001;
N = 3;
true = 1;
false = 0;
type myarray = array [1..N] of real;
var x, xE: real;
y, yE, z: myarray;
function FUNC( which: integer; x, y: real ): real;
begin
case which of
1: FUNC := 2 * x;
2: FUNC := 3 * x * x;
3: FUNC := y;
otherwise
FUNC := INVALID_NUM
end {case}
end {FUNC};
procedure FKT( x: real; y: myarray; n: integer; var z: myarray );
var i: integer;
begin
for i:= 1 to n do z[i] := FUNC( i, x, y[i] )
end {FKT};
procedure RK1ST( x: real; y: myarray; h: real; n: integer;
var xe: real; var ye: myarray );
var j, k: integer;
w: myarray;
a: array [1..5] of real;
begin
a[1] := h/2; a[2] := h/2; a[5] := h/2; a[3] := h; a[4] := h;
xe := x;
for k:= 1 to n do
begin
ye[k] := y[k];
w[k] := y[k]
end {for k};
for j:= 1 to 5 do
begin
FKT( xe, w, n, z );
xe := x + a[j];
for k:= 1 to n do
begin
w[k] := y[k] + a[j] * z[k];
ye[k] := ye[k] + a[j+1] * z[k] / 3
end {for k}
end {for j}
end {RK1ST};
procedure RK( x: real; y: myarray; n: integer; eps: real;
xE: real; var yE: myarray; fi: integer );
var y1, y2, y3: myarray;
x1, x2, x3, H: real;
s, Hs: real;
out, j, k: integer;
do_BB, do_EE: integer; { dummy flags }
begin
if fi=true then
begin
H := xE - x;
s := 0
end {if}
else
begin
H := Hs;
out := false
end {else};
while out=false do
begin
if ((x + 2.01 * H - xE) > 0) = (H > 0) then
begin
Hs := H;
out := true;
H := ( xE - x ) / 2
end {if};
RK1ST( x, y, 2 * H, n, x1, y1 );
do_BB := true;
{ BB }
while do_BB=true do
begin
do_EE := true;
RK1ST( x, y, H, n, x2, y2 );
RK1ST( x2, y2, H, n, x3, y3 );
for k:= 1 to n do
begin
if ( y1[k] - y3[k] ) > eps then
begin
{ CC }
H := 0.5 * H;
out := false;
x1 := x2;
for j:= 1 to n do y1[j] := y2[j];
do_EE := false;
end {if}
end {for k};
{ EE - second part of BB }
if do_EE=true then
begin
x:= x3;
if out=true then
for k:= 1 to n do yE[k] := y3[k]
else
begin
for k:= 1 to n do y[k] := y3[k];
if s = 5 then
begin
s := 0;
H := 2 * H
end {if s};
s := s + 1;
do_BB := false;
end {else}
end {if do_EE}
end {while do_BB};
end {while out};
end {RK};
begin
RK( x, y, N, EPSILON, xE, yE, true )
end.
--------------------------------------------------
Discussion of the Runge-Kutta Algorithm
written in Algol60, C and Pascal
(1) Algol60 --> C
* Since C doesn't support function/procedure
nesting, we have to change the nested Algol60
function definitions to functions defined
in parallel.
Also, because of this reason, I needed to add
one more parameter in procedure RK1ST(.., int n, ..)
in order to make this function more general instead
of using the global defined "N".
* The function parameters in Algol60 are passed
by default "call-by-name" or by "call-by-value"
if specified. In C, we need to use "call-by-value"
and "call-by-reference". For "call-by-reference"
in C, we need to use pointer to pass the parameters.
* Algol60 uses "goto" statements. Although it is easy
to follow the flow, sometimes they are confusing.
In C, I changed those goto's to several "while" and
"if" statements.
* Compared to Algol60, the code written in C is much
cleaner.
* For procedure RK(..), some parameters used in
Algol60 version are not necessary. So, my change
here is:
Algol60: RK(x,y,n,FKT,eps,eta,xE,yE,fi) -->
C: RK(x,y,n,eps,xE,yE,fi)
(2) Algol60(or C) --> Pascal
* Pascal also supports function/procedure nesting.
I didn't use nesting for the conversion.
* Compared to C, I feel more comfortable with the
parameter passing in Pascal than in C. In C, we
have to use pointer for "call-by-reference" while
in Pascal we use "var" to specify that.
* Compared to C again, C has "break" statement
which sometimes might be very efficient to get
out from the loops. Pascal doesn't have this
property. At one place in my C code, I used a
"break" statement.
* As in the C version, I also got rid of all the
"goto"s, and I used a lot of "while"s and "if"s.
* For the function/procedure parameters, I did the
same changes here as what I did for the C version.