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