/********************************************************** * SOLVING DIFFERENTIAL EQUATIONS OF ORDER 1 * * (Prediction-correction method) * * ------------------------------------------------------- * * Reference: * * "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc * * DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991" * * [BIBLI 03]. * * * * C++ Version By J-P Moreau, Paris * * (www.jpmoreau.fr) * * ------------------------------------------------------- * * SAMPLE RUN: * * (Integrate y'=(4xy+4x*sqrt(y))/(1+x²) from x=0 to x=10 * * with initial condition y(0)=1) * * * * Input x begin: 0 * * Input x end : 10 * * Input y begin: 1 * * Input number of points: 10 * * Input finesse: 10 * * * * X Y * * ----------------------------- * * 1.000000 8.999984 * * 2.000000 80.999881 * * 3.000000 360.999496 * * 4.000000 1088.998512 * * 5.000000 2600.996484 * * 6.000000 5328.992838 * * 7.000000 9800.986875 * * 8.000000 16640.977767 * * 9.000000 26568.964559 * * 10.000000 40400.946171 * * ----------------------------- * * * **********************************************************/ #include #include #define MAXDATA 100 double h, tx[MAXDATA],ty[MAXDATA], x,xi,xf,yi; int fi,i,n; // y'=(4xy+4x*sqrt(y))/(1+x²) double f(double x, double y) { return (4*x*y+4*x*sqrt(y))/(1+x*x); } //classical Runge-Kutta method of order 4 void runge_kutta(double *y) { double a,b,c,d; a=h*f(x,*y); b=h*f(x+h/2,*y+a/2); c=h*f(x+h/2,*y+b/2); x+=h; d=h*f(x,*y+c); *y =*y + (a+b+b+c+c+d)/6; } /********************************************************* * Prediction-correction method * * ------------------------------------------------------ * * INPUTS: * * xi begin x value * * xf end x value * * y1 begin y value (at x=xi) * * m number of points to calculate * * fi finesse (number of intermediate * * points (for example 10) * * OUTPUTS: * * tx table of x values (m values) * * ty table of y values (m values) * * * * DESCRIPTION: * * The algorithm has the following steps: * * 1. calculate y2, y3, y4 using a classical Runge- * * Kutta method of order 4. * * 2. from point (x4, y4), first estimate y(n+1) by * * formula: * * y(n+1)=y(n) + H/24(55y'(n)-59y'(n-1)+37y'(n-2) * * -9y'(n-3) * * then continue with formula: * * y(n+1)=y(n) + H/24(9y'(n+1)+19y'(n)-5y'(n-1) * * +y'(n-2), * * noting that y'(n+1)=f(x(n+1),y(n+1)) with the * * estimated value of y(n+1), until convergence is * * obtained. * *********************************************************/ void equadiff_pc(double *tx,double *ty,double xi,double xf, double yi,int m,int fi) { double z=yi,y,w,p[4]; int i,j,k; long ni; if ((m>MAXDATA) || (fi<1)) return; h = (xf-xi)/fi/m; p[3]=f(xi,yi); tx[0]=xi; ty[0]=yi; for (k=0, i=1; i<=m; i++) { ni=(long) (i-1)*fi-1; for (j=1; j<=fi; j++) { x=xi+h*(ni+j); if (++k<4) { runge_kutta(&z); p[3-k]=f(x,z); } else { x+=h; w=z+h/24*(55*p[0]-59*p[1]+37*p[2]-9*p[3]); do { y=w; w=z+h/24*(9*f(x,y)+19*p[0]-5*p[1]+p[2]); } while (fabs(y-w) > 1e-10); z=w; p[3]=p[2]; p[2]=p[1]; p[1]=p[0]; p[0]=f(x,z); } } // j loop tx[i]=x; ty[i]=z; } // i loop } void main() { printf("\n Input x begin: "); scanf("%lf",&xi); printf("\n Input x end : "); scanf("%lf",&xf); printf("\n Input y begin: "); scanf("%lf",&yi); printf("\n Input number of points: "); scanf("%d",&n); printf("\n Input finesse: "); scanf("%d",&fi); equadiff_pc(tx,ty,xi,xf,yi,n,fi); //print results printf("\n"); printf(" X Y \n"); printf(" -----------------------------\n"); for (i=1; i<=n; i++) printf("%12.6f %12.6f\n",tx[i],ty[i]); printf(" -----------------------------\n"); } //end of file teqdifpc.cpp