#include
#include
/* 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));
}
}