0
/*
---
Algo9-8.c C program for implementing Algorithm 9.8
Algorithm translated to C by: Dr. Norman Fahrer
IBM and Macintosh verification by: Daniel Mathews
NUMERICAL METHODS: C Programs, (c) John H. Mathews 1995
To accompany the text:
NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions: ISBN 0-13-625047-5
This free software is compliments of the author.
E-mail address: in%"mathews@fullerton.edu"
Algorithm 9.8 (The Hamming Method).
Section 9.6, Predictor-Corrector Method, Page 473
---
•
/
/*
---
Algorithm 9.8 (The Hamming Method)
To approximate the solution of the initial value problem y' = f(t,y)
with y(a) = y_0 over [a,b] by using the predictor
4h
p_(k+1) = y_(k-3) + --- [2*f_(k-2) - f_(k-1) + 2*f_k]
3
and the corrector :
-y_(k-2) + 9 y_k 3h
y_(k+1) = --- + --- [-f_(k-1) + 2*f_k + f_(k+1)].
8 3
User has to supply functions named : ffunction
An example is included in this program.
---
•
/
1.
include<stdio.h>
2.
include<stdlib.h>
3.
include<math.h>
4.
define MAX 500
/* define prototype for USER-SUPPLIED function f(x) */
double ffunction(double t, double y);
/* EXAMPLE for "ffunction" */
double ffunction(double t, double y)
{
return ( (t - y) / 2.0 );
}
/* --- */
/* Main program for algorithm 9.8 */
void main(void)
{
int I, K; /* Loop counter */
int N; /* Number of steps > 3 */
double A, B, Y[MAX]; /* Endpoints and inital value */
double H; /* Compute the step size */
double T[MAX];
double F1, F2, F3, F4; /* Function values */
double Hmin, Hmax; /* Minimum and maximum step size */
double Pold, Cold, Cnew; /* Predictor and Corrector */
double Pmod, Pnew; /* Modifier a Predictor */
printf("---The Hamming Method---n");
printf("--- Example 9.13 on page 468 ---n");
printf("---n");
printf("Please enter endpoints of the interval [A,B]:n");
printf("For used Example type : 0, 3.0n");
scanf("%lf, %lf", &A, &B);
printf("You entered [%lf, %lf]n", A, B);
printf("---n");
printf("Please enter number of steps: ( > 3 and < %d !)n", MAX+1);
scanf("%d",&N);
if( (N > MAX) || (N < 4) )
{
printf(" Number of steps must be greater than 3 and less than %dn",MAX+1);
printf(" Terminating. Sorryn");
exit(0);
}
printf("---n");
printf("You need all together FOUR initial values. You cann");
printf("use the Runge-Kutta method to compute the 2nd, 3rd andn");
printf("4th from the 1st one.n");
printf("Please enter initial values Y[0], Y[1], Y[2], Y[3] :n");
printf("Example 9.13 page 468: 1, 0.94323919, 0.89749071, 0.86208736n");
scanf("%lf, %lf, %lf, %lf", &Y[0], &Y[1], &Y[2], &Y[3]);
printf("You entered Y[0] = %lfn", Y[0]);
printf("You entered Y[1] = %lfn", Y[1]);
printf("You entered Y[2] = %lfn", Y[2]);
printf("You entered Y[3] = %lfn", Y[3]);
/* Compute the step size and initialize */
H = (B - A) / N;
T[0] = A;
for( K = 1; K <= 3; K++) T[K] = A + K * H;
F1 = ffunction(T[1],Y[1]);
F2 = ffunction(T[2],Y[2]);
F3 = ffunction(T[3],Y[3]);
Pold = 0;
Cold = 0;
for( K = 3; K <= N-1; K++)
{
/* Milne Predictor */
Pnew = Y[K-3] + 4.0 * H * (2.0*F1 - F2 + 2.0*F3) / 3.0;
/* Modifier */
Pmod = Pnew + 112.0 * (Cold - Pold) / 121.0;
/* Next mesh point */
T[K+1] = A + H * (K+1);
/* Evaluate f(t,y) */
F4 = ffunction(T[K+1],Pmod);
/* Hamming Corrector */
Cnew = (9.0 * Y[K] - Y[K-2] + 3.0 * H * ( -F2 + 2.0 * F3 + F4))/8.0;
/* New value at Y_(k+1) */
Y[K+1] = Cnew + 9.0 * (Pnew - Cnew) / 121.0;
/* Update the values */
Pold = Pnew;
Cold = Cnew;
F1 = F2;
F2 = F3;
F3 = ffunction(T[K+1], Y[K+1]);
}
/* Output */
for ( K = 0; K <= N; K++)
{
printf("K = %d, T[K] = %lf, Y[K] = %lfn", K, T[K], Y[K]);
}
} /* End of main program */
Tümünü Göster
