! ! =====> 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*