/************************************************************ * Differential equation y"=f(x,y,y') by Stormer's method * * --------------------------------------------------------- * * SAMPLE RUN: * * Integrate y" = 8yy / (1 + 2x) from x=0 to x=1, * * with initial conditions: x(0)=0, y(0)=1 and y'(0)=-2 * * and compare with exact solution: y = 1 / (1 + 2x) * * * * Output file (stormer.lst): * * * * --------------------------------------------------------- * * Differential equation y"=f(x,y,y') by Stormer's method * * --------------------------------------------------------- * * X Y Y exact Error * * 0.000 1.000000 1.000000 0.0000000000 * * 0.010 0.980392 0.980392 0.0000000001 * * 0.020 0.961538 0.961538 0.0000000295 * * 0.030 0.943396 0.943396 0.0000000457 * * 0.040 0.925926 0.925926 0.0000000974 * * 0.050 0.909091 0.909091 0.0000001285 * * ... ... ... ... * * 0.950 0.344866 0.344828 0.0000381695 * * 0.960 0.342505 0.342466 0.0000388874 * * 0.970 0.340176 0.340136 0.0000396196 * * 0.980 0.337878 0.337838 0.0000403406 * * 0.990 0.335612 0.335570 0.0000410721 * * 1.000 0.333375 0.333333 0.0000418231 * * * * End of file. * * C++ Version By J-P Moreau * * (www.jpmoreau.fr) * ************************************************************/ #include #include // Example: y" = 8 yy / (1 + 2x) double F(double x,double y,double z) { return z; } double G(double x,double y,double z) { return (8.0*y*y/(1.0+2.0*x)); } // exact solution double Fx(double x) { return (1.0/(1.0+2.0*x)); } void RK4D2(double x,double y,double z,double h, double *x1,double *y1,double *z1) { double c1,d1,c2,d2,c3,d3,c4,d4; c1=F(x,y,z); d1=G(x,y,z); c2=F(x+h/2,y+h/2*c1,z+h/2*d1); d2=G(x+h/2,y+h/2*c1,z+h/2*d1); c3=F(x+h/2,y+h/2*c2,z+h/2*d2); d3=G(x+h/2,y+h/2*c2,z+h/2*d2); c4=F(x+h,y+h*c3,z+h*d3); d4=G(x+h,y+h*c3,z+h*d3); *x1=x+h; *y1=y+h*(c1+2*c2+2*c3+c4)/6; *z1=z+h*(d1+2*d2+2*d3+d4)/6; } void main() { double c[5],x[5],y[5],z[5]; double h, a1, a2, a3; double yex,er; int k; FILE *fp; h=0.01; a1=1.08333333333333; a2=-2.0*(a1-1.0); a3=a1-1.0; fp=fopen("stormer.lst","w"); fprintf(fp,"-----------------------------------------------------------------\n"); fprintf(fp," Differential equation y\"=f(x,y,y') by Stormer's method \n"); fprintf(fp,"-----------------------------------------------------------------\n"); // initial conditions x[1]=0.0; y[1]=1.0; z[1]=-2.0; yex=Fx(x[1]); er=0.0; fprintf(fp," X Y Y exact Error\n"); fprintf(fp," %f %f %f %f\n",x[1],y[1],yex,er); for (k=1; k<3; k++) { //call Runge-Kutta for first 2 steps RK4D2(x[k],y[k],z[k],h,&x[k+1],&y[k+1],&z[k+1]); yex=Fx(x[k+1]); er=fabs(yex-y[k+1]); fprintf(fp," %f %f %f %f\n",x[k+1],y[k+1],yex,er); } //main Stormer loop while (x[3]<1.0) { // end x value = 1.0 for (k=2; k<5; k++) { c[k]=G(x[5-k],y[5-k],z[5-k]); } y[4]=2.0*y[3]-y[2]+h*h*(a1*c[2]+a2*c[3]+a3*c[4]); x[4]=x[3]+h; yex=Fx(x[4]); er=fabs(yex-y[4]); fprintf(fp," %f %f %f %f\n",x[4],y[4],yex,er); for (k=1; k<4; k++) { x[k]=x[k+1]; y[k]=y[k+1]; } } fprintf(fp,"\n END OF FILE."); fclose(fp); printf("\n Results in file stormer.lst\n\n"); } //end of file stormer.cpp