/************************************************************************* * Solve a boundary value problem for a second order DE of the form: * * y" - p(x)*y = f(x) with conditions: y(x0)=a and y(x1)=b (x0<=x<=x1) * * using a Runge-Kutta integration method. * * ---------------------------------------------------------------------- * * Explanations: if y0(x) is a particular solution of y"-p(x)y = f(x) (1) * * and y1(x), y2(x) two lineary independant solutions of y"-p(x)y = 0 (2) * * (right hand side removed), then the general solution of (1) has the * * form: y(x) = y0(x) + c1*y1(x) + c2*y2(x), the constants c1, c2 are * * obtained from boundary conditions. * * A more efficient way to proceed is: 1) calculate a particular solution * * y1(x) of (2) such as y1(x0)=0 and 2) calculate a particular solution * * y2(x) of (1) such as y2(x0)=a, then the general solution of (1) has the* * form: y(x) = y1(x) + C * y2(x), constant C is determined by: * * C * y2(x1) + y1(x1) = b * * ---------------------------------------------------------------------- * * SAMPLE RUN: * * Integrate y" - (1/x)y' - (3/x²)y = 5x² from x0=1 to x1=2 with boundary * * conditions: y(x0)=11 and y(x1)=36. The exact solution is: * * y(x) = 8/x + 2x² + x4 * * The initial 2nd order problem is replaced by the following first order * * DE systems: * * 1st system 2nd system * * y' = z y' = z * * z' = z/x² + 3y/x² z' = z/x² + 3y/x² + 5x² * * with z'(x0)=0 with z'(x0)=11 * * * * X Y estimated Y true Absolute Error * * -------------------------------------------------- * * 1.0000 11.0000000 11.0000000 0.00e+000 * * 1.0500 11.1498044 11.1498039 5.73e-007 * * 1.1000 11.3988283 11.3988273 9.96e-007 * * 1.1500 11.7472793 11.7472780 1.31e-006 * * 1.2000 12.1962682 12.1962667 1.53e-006 * * 1.2500 12.7476579 12.7476563 1.69e-006 * * 1.3000 13.4039479 13.4039462 1.78e-006 * * 1.3500 14.1681840 14.1681822 1.84e-006 * * 1.4000 15.0438876 15.0438857 1.85e-006 * * 1.4500 16.0349995 16.0349976 1.8E-0006 * * 1.5000 17.1458351 17.1458333 1.8E-0006 * * 1.5500 18.3810483 18.3810466 1.7E-0006 * * 1.6000 19.7456016 19.7456000 1.6E-0006 * * 1.6500 21.2447426 21.2447411 1.5E-0006 * * 1.7000 22.8839837 22.8839824 1.3E-0006 * * 1.7500 24.6690860 24.6690848 1.1E-0006 * * 1.8000 26.6060454 26.6060444 9.6E-0007 * * 1.8500 28.7010813 28.7010806 7.5E-0007 * * 1.9000 30.9606268 30.9606263 5.2E-0007 * * 1.9500 33.3913206 33.3913204 2.7E-0007 * * 2.0000 36.0000000 36.0000000 5.0E-0014 * * * * ---------------------------------------------------------------------- * * REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en Basic * * et en Pascal By Claude Nowakowski, Edition du P.S.I., 1984"* * [BIBLI 04]. * * * * C++ Release By J-P Moreau, Paris. * *************************************************************************/ #include #include #define SIZE 100 double X[SIZE]; double Y[3][SIZE], Z[3][SIZE]; double cc,er,xl,x0,h,ya,yb,yex; int k,kl,l; // y'=z double F(double x, double y,double z) { return z; } // l=1: z'=z/x²+3y/x² or l=2: z'=z/x²+3y/x²+5x² double G(double x,double y,double z) { double gg; gg=z/x+3.0*y/(x*x); if (l==2) gg += 5.0*x*x; return gg; } // exact solution: y=8/x + 2x3 + x4 double FX(double x) { return 8.0/x+x*x*x*(x+2.0); } // Integrate sytem from x to x+h using Runge-Kutta void RK4(double x,double y,double z,double h,double *x1,double *y1, double *z1) { double c1,c2,c3,c4,d1,d2,d3,d4,h2; c1=F(x,y,z); d1=G(x,y,z); h2=h/2.0; c2=F(x+h2,y+h2*c1,z+h2*d1); d2=G(x+h2,y+h2*c1,z+h2*d1); c3=F(x+h2,y+h2*c2,z+h2*d2); d3=G(x+h2,y+h2*c2,z+h2*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.0*c2+2.0*c3+c4)/6.0; *z1=z+h*(d1+2.0*d2+2.0*d3+d4)/6.0; } void main() { x0=1.0; //starting x xl=2.0; //ending x kl=20; //number of steps in x h=(xl-x0)/kl; //integration step ya=11.0; //starting value for y yb=36.0; //final value for y'' X[0]=x0; // integration loop for (l=1; l<3; l++) { Y[l][0]=(l-1)*ya; Z[l][0]=1.0; for (k=0; k