/******************************************************************* * Differential equations of order N * * by Runge-Kutta method of order 4 * * with boundary conditions * * ---------------------------------------------------------------- * * Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0" by Marc * * DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991. * * * * C++ version By J-P Moreau, Paris * * (www.jpmoreau.fr) * * ---------------------------------------------------------------- * * SAMPLE RUN: * * * * Example: integrate y"=(2y"y"+yy)/y from x=4 to x=6 * * with initial conditions: y(4)=2 and y'(4)=-2 * * and final condition y(6)=2 * * * * Exact solution is: y = 2cos(1)/cos(x-5) * * * * DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER N * * of type y(n) = f(y(n-1),...,y"",y,x) * * with boundary conditions * * * * order of equation: 2 * * * * begin value x : 4 * * end value x : 6 * * y0 value at x0 : 2 * * y1 value at x0 : -2 * * Desired value at xf: 2 * * number of points : 11 * * finesse : 20 * * Max. iterations : 100 * * * * Last iteration: * * u =-3.114816 du =-0.000001 * * dyb=0.000000 67 steps * * * * X Y * * --------------------------- * * 4.000000 2.000000 * * 4.200000 1.551018 * * 4.400000 1.309291 * * 4.600000 1.173217 * * 4.800000 1.102583 * * 5.000000 1.080605 * * 5.200000 1.102583 * * 5.400000 1.173217 * * 5.600000 1.309291 * * 5.800000 1.551018 * * 6.000000 2.000000 * * * *******************************************************************/ #include #include #define FALSE 0 #define TRUE 1 int fi,i,ndata,ordre; double xi,xf; double yi[10]; char s[2]; double du,dyb,olddyb,u,yb; int maxiter, niter; // Example: y"=(2y'y'+yy)/y double fp(double x, double *y) { if (fabs(y[0])<1e-12) y[0]=1e-12; return (2*y[1]*y[1]+y[0]*y[0])/y[0]; } void Display(double *t,int m,double xi,double xf) { double h,x; int i; h=(xf-xi)/(m-1); x=xi-h; printf("\n X Y \n"); printf("------------------------\n"); for (i=1; i1e-7 && niter fabs(olddyb)) // overshoot? if (niter==1) du =-du; else du =-du/2.0; olddyb =dyb; niter++; } printf("\n Last iteration:\n"); printf(" u =%f du =%f\n", u, du); printf(" dyb=%f %d steps\n", dyb, niter); // print final results equadiffn(vect,xi,xf,yi,ndata,ordre,fi,TRUE); printf("\n Any key + Enter to end program..."); scanf("%s",s); } // End of file tshoot.cpp