!
! =====> Program - P123.F90
!
PROGRAM DifDmo
IMPLICIT NONE
! Program to demonstrate numerical solutions
! to ordinary, first order differential equations.
INTEGER :: I, N
REAL :: X, Xf, h, &
Y0, Y_Eul, Y_RK4
PRINT *, 'This is Program >> P123 = Differential Equations'
!
!
PRINT *, 'First order ODEs F = -y*y for example'
PRINT *, 'Enter final (high) x value: '
READ *, Xf
Print *, Xf
PRINT *, 'Number of steps to take: '
READ *, N
Print *, N
! Step size, initial X, Y (given)
h = Xf/N
X = 0.0
Y0 = 1
PRINT *, 'Step size: ', h
PRINT 101, 'Step', 'x value', 'Euler', 'RK-4'
! Compare Euler and Runge-Kutta 4'th order
! at each step.
Y_Eul = Y0
Y_RK4 = Y0
L1: DO I = 1, N
Y_Eul = Euler(X, Y_Eul, h)
Y_RK4 = RngKt4(X, Y_RK4, h)
X = X + h
PRINT 103, I, X, Y_Eul, Y_RK4
END DO L1
STOP
101 FORMAT(//' ', A5, 6A12)
103 FORMAT(' ', I5, 3F12.8)
END PROGRAM DifDmo
FUNCTION F(x, y)
! The function to be integrated.
! I.e. y' = dy/dx = F(x, y)
F = x ! Just to keep FTN90 quiet
F = -y*y
RETURN
END FUNCTION F
! ALGORITHMS
FUNCTION Euler(x, y, h)
! Uses F(x, y) to advance the solution
! from x to x+h.
REAL :: x, y, h
Euler = y + h*F(x, y)
RETURN
END FUNCTION Euler
FUNCTION RngKt4(x, y, h)
! Uses F(x, y) to advance the solution
! from x to x+h.
REAL :: x, y, h, &
h_hf, &
x_hf, &
y_hf_1, &
y_hf_2, &
y_1
h_hf = h/2
x_hf = x + h_hf
! Form various estimates of the solution
! at x+h/2 and x+h, then combine them
! using the formula below.
y_hf_1 = Euler(x, y, h_hf)
y_hf_2 = Euler(x_hf, y_hf_1, h_hf)
y_1 = Euler(x_hf, y_hf_2, h)
RngKt4 = y + h/3 * &
((f(x, y) + f(x+h, y_1))/2 + &
f(x_hf, y_hf_1) + f(x_hf, y_hf_2))
RETURN
END FUNCTION RngKt4
#include <stdio.h>
#include <math.h>
/* DIFFDEMO.C */
/* Program to demonstrate numerical solutions to ordinary, first order
differential equations (initial value problems). Starting from
(x, y) = (0, y0), the solution to y' = dy/dx is advanced through
N steps of size h to a final point (Xf, Yf).
Two "self-starting" methods are used here. The primitive Euler method
introduces the topic, while the 4'th order Runge-Kutta formula is
useful in practice.
Note that many other more sophisticated and accurate methods of numerical
solution of differential equations exist.
*/
double y_solution(double x, double y0){
/* The analytical solution y = y(x) to the differential equation
y' = F(x, y) with initial condition y(0) = y0. */
double solution;
solution = 1/(x + 1/y0);
return(solution);
}
double f(double x, double y){
/* The function to be integrated. i.e. y' = dy/dx = F(x, y) */
double f;
f= - y*y;
return(f);
}
double euler(double x, double y, double h){
/* Uses F(x, y) to advance the solution of the first order
ordinary differential equation from x to x+h. Note that
the "truncation error" in the first order Euler method is
O(h), so the estimated solution is not accurate. */
return( y+ h*f(x,y) );
}
double runge_kutta4(double x, double y, double h){
/* Uses F(x, y) to advance the solution of the first order
ordinary differential equation from x to x+h. The truncation
error is O(h^4), so the method is reasonable accurate despite
the simplicity of coding. */
double t1, t2, t3, t4;
t1= h * f(x,y);
t2= h * f(x+h/2.0, y+t1/2.0);
t3= h * f(x+h/2.0, y+t2/2.0);
t4= h * f(x+h, y+t3);
y= y + (t1 + 2.0*t2 + 2.0*t3 + t4) /6.0;
return(y);
}
main(){
int a,i,n;
double x, xf, h, y0, y_true, y_euler, y_rk4;
printf("Numerical Solutions to differential equations\n\n\n");
printf(" This program will use the Euler method as well as the much\n");
printf("better Runge-Kutta order 4 method to approximate the solution.\n");
printf("The diferential equation used is y' = - (y-4) / x, with the \n");
printf("initial condition being y(1) = 0. You will be prompted for \n");
printf("the final value of x, and the number of steps. The two \n");
printf("methods will be compared at each step with the analytical \n");
printf("solution, so that you can compare their accuracy.\n\n\n");
do{
printf("Enter final (high) x value: ");
scanf("%lf",&xf);
} while(xf < 0);
do{
printf("Enter number of steps to take: ");
scanf("%d",&n);
} while(n < 0);
x= 0.0; /* Initial X */
h= (xf-x)/n; /* Step size */
y0= 1.0; /* Initial Y (given as part of initial value problem) */
printf("\nThe differential equation will be solved\n");
printf("from x = 1 to x = %7.2lf\n",xf);
printf("using %5d steps of size %7.5lf\n\n\n",n,h);
printf(" Step x-value true-y Euler error RK-4 error\n");
for(a=0;a<=60;a++)
printf("-");
printf("\n");
y_euler = y0;
y_rk4 = y0;
/* Compare true solution with Euler and Runge-Kutta 4'th order
at each step */
for(i=0;i<=n;i++){
y_euler= euler(x, y_euler, h);
y_rk4= runge_kutta4(x, y_rk4, h);
x= x + h;
y_true= y_solution(x, y0);
printf("%5d %7.4lf %7.4lf %7.4lf %7.4lf %7.4lf %7.4lf\n",
i, x, y_true, y_euler, (y_euler - y_true), y_rk4, (y_rk4 - y_true));
}
}
PROGRAM DiffDemo;
{ Program to demonstrate numerical solutions to ordinary, first order
differential equations (initial value problems). Starting from
(x, y) = (0, y0), the solution to y' = dy/dx is advanced through
N steps of size h to a final point (Xf, Yf).
Two "self-starting" methods are used here. The primitive Euler method
introduces the topic, while the 4'th order Runge-Kutta formula is
useful in practice.
Note that many other more sophisticated and accurate methods of numerical
solution of differential equations exist.
This program was compiled with Turbo Pascal V5.5, but should be
compilable with Turbo Pascal versions 3.0 and above. Deviations from
standard Pascal are also minimal, particularly in the coding of the
algorithms themselves. }
VAR
I, N: Integer;
X, Xf, h,
Y0, Y_True,
Y_Euler,
Y_RK4: Real;
FUNCTION F(x, y: Real): Real;
{ The function to be integrated. I.e. y' = dy/dx = F(x, y) }
BEGIN
F:= y;
END; { FUNCTION F }
FUNCTION Y_solution(x, y0: Real): Real;
{ The analytical solution y = y(x) to the differential equation
y' = F(x, y) with initial condition y(0) = y0. }
BEGIN
Y_solution:= exp(x) + Y0;
END; { FUNCTION Y_solution }
(* FUNCTION F(x, y: Real): Real;
{ The function to be integrated. I.e. y' = dy/dx = F(x, y) }
BEGIN
F:= -y*y;
END; { FUNCTION F }
FUNCTION Y_solution(x, y0: Real): Real;
{ The analytical solution y = y(x) to the differential equation
y' = F(x, y) with initial condition y(0) = y0. }
BEGIN
Y_solution:= 1/(x + 1/y0);
END; { FUNCTION Y_solution } *)
FUNCTION Euler(x, y, h: Real): Real;
{ Uses F(x, y) to advance the solution of the first order
ordinary differential equation from x to x+h. Note that
the "truncation error" in the first order Euler method is
O(h), so the estimated solution is not accurate. }
BEGIN
Euler:= y + h*F(x, y);
END; { FUNCTION Euler }
FUNCTION Runge_Kutta4(x, y, h: Real): Real;
{ Uses F(x, y) to advance the solution of the first order
ordinary differential equation from x to x+h. The truncation
error is O(h^4), so the method is reasonable accurate despite
the simplicity of coding. }
VAR
h_half,
x_half,
y_half_1,
y_half_2,
y_1: Real;
BEGIN
h_half := h/2;
x_half := x + h_half;
{ Form various estimates of the solution at x+h/2 and x+h, then
combine them using the formula below. }
y_half_1 := Euler(x, y, h_half);
y_half_2 := Euler(x_half, y_half_1, h_half);
y_1 := Euler(x_half, y_half_2, h);
Runge_Kutta4:= y + h/3 *
((f(x, y) + f(x+h, y_1))/2 +
f(x_half, y_half_1) + f(x_half, y_half_2));
END; { FUNCTION Runge_Kutta4 }
BEGIN { Main }
{ Get input from user---target x, Xf, and number of steps, N }
REPEAT
Write('Enter final (high) x value: ');
ReadLn(Xf);
UNTIL (Xf > 0);
REPEAT
Write('Enter number of steps to take: ');
ReadLn(N);
UNTIL (N > 0);
h:= Xf/N; { Step size }
X:= 0.0; { Initial X }
Y0:= 1; { Initial Y (given as part of initial value problem) }
WriteLn; WriteLn;
WriteLn('The differential equation will be solved ');
WriteLn('from x = 0 to x = ', Xf:7:2);
WriteLn('using ', N:1, ' steps of size ', h:7:4);
WriteLn; WriteLn;
WriteLn('Step':5, 'x value':10, 'True y':10,
'Euler':10, 'Error':10, 'RK-4':10, 'Error':10);
Write(0:5, X:10:4, Y0:10:4);
FOR I:= 1 TO 4 DO Write('--':10); WriteLn;
Y_Euler:= Y0;
Y_RK4:= Y0;
{ Compare true solution with Euler and Runge-Kutta 4'th order
at each step }
FOR I:= 1 TO N DO
BEGIN
Y_Euler:= Euler(X, Y_Euler, h);
Y_RK4:= Runge_Kutta4(X, Y_RK4, h);
X:= X + h;
Y_True:= Y_solution(X, Y0);
WriteLn(I:5, X:10:4, Y_True:10:4,
Y_Euler:10:4, (Y_Euler - Y_True):10:4,
Y_RK4:10:4, (Y_RK4 - Y_True):10:4);
END; { FOR I }
END.
Last modified: 08/07/97