/******************************************************************* * Differential equations with p variables of order 1 * * by Runge-Kutta method of order 4 * * ---------------------------------------------------------------- * * Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc * * DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991" * * [BIBLI 03]. * * * * C++ Version 1.1 By J-P Moreau, Paris * * (www.jpmoreau.fr) * * ---------------------------------------------------------------- * * SAMPLE RUN: * * * * Example #1: integrate system of equations from x=0 to x=3: * * y1' = y2 + y3 - 3*y1 * * y2' = y1 + y3 - 3*y2 * * y3' = y1 + y2 - 3*y3 * * with initial conditions: y1(0)=1, y2(0)=2 and y3(0)=-1 * * * * * * DIFFERENTIAL EQUATION WITH P VARIABLE OF ORDER 1 * * of type yi' = f(y1,y2,...,yn), i=1..n * * * * number of variables: 3 * * begin value x : 0 * * end value x : 3 * * y1 value at x0 : 1 * * y2 value at x0 : 2 * * y3 value at x0 : -1 * * number of points : 7 * * finesse : 30 * * * * X Y1 Y2 Y3 * * -------------------------------------------------------- * * 0.000000 1.000000 2.000000 -1.000000 * * 0.500000 0.449466 0.584801 0.178795 * * 1.000000 0.251358 0.269674 0.214727 * * 1.500000 0.149580 0.152058 0.144622 * * 2.000000 0.090335 0.090671 0.089664 * * 2.500000 0.054738 0.054784 0.054648 * * 3.000000 0.033193 0.033200 0.033181 * * -------------------------------------------------------- * * * * Example #2: integrate system of equations from x=0 to PI*Sqrt(2) * * y1' = y2 * * y2' = -4*y1 - 3*y3 * * y3' = y4 * * y4' = -8*y1 - 2*y3 * * with initial conditions: y1(0)=3, y2(0)=0, y3(0)=4, y4(0)=0 * * * * * * DIFFERENTIAL EQUATION WITH P VARIABLE OF ORDER 1 * * of type yi' = f(y1,y2,...,yn), i=1..n * * * * number of variables: 4 * * begin value x : 0 * * end value x : 4.442883 * * y1 value at x0 : 3 * * y2 value at x0 : 0 * * y3 value at x0 : 4 * * y4 value at x0 : 0 * * number of points : 9 * * finesse : 30 * * * * X Y1 Y2 Y3 Y4 * * -------------------------------------------------------- * * 0.000000 3.000000 0.000000 4.000000 0.000000 * * 0.555360 0.000000 -8.485281 0.000000 -11.313708 * * 1.110721 -3.000000 -0.000001 -4.000000 -0.000002 * * 1.666081 -0.000001 -8.485281 -0.000001 11.313708 * * 2.221442 3.000000 0.000003 4.000000 0.000003 * * 2.776802 0.000001 -8.485281 0.000002 -11.313708 * * 3.332162 -3.000000 -0.000004 -4.000000 -0.000005 * * 3.887523 -0.000002 8.485281 -0.000002 11.313708 * * 4.442883 3.000000 0.000005 4.000000 0.000007 * * -------------------------------------------------------- * * * * Release 1.1: Added example #2 (Sept. 2008). * *******************************************************************/ #include #include #define SIZE 50 int fi,i,ndata,p; double xi,xf; double yi[10]; double v1[SIZE],v2[SIZE],v3[SIZE],v4[SIZE]; //Example: y1'=y2+y3-3y1, y2'=y1+y3-3y2, y3'=y1+y2-3y3 double fp(int k, double x, double *y) { switch(k) { case 0: return y[1]+y[2]-3*y[0]; case 1: return y[0]+y[2]-3*y[1]; case 2: return y[0]+y[1]-3*y[2]; default: return 0; } } /*Exemple : y1'=y2, y2'=-4y1-3y3, y'3=y4, y'4=-8y1-2y3 double fp(int k, double x, double *y) { switch(k) { case 0: return y[1]; case 1: return (-4*y[0]-3*y[2]); case 2: return y[3]; case 3: return (-8*y[0]-2*y[2]); default: return 0; } } */ // print results (p>1) void Display(double *t1,double *t2,double *t3,double *t4,int m,int p,double xi,double xf) { int i; double h,x; h=(xf-xi)/(m-1); x=xi-h; printf("\n X"); for (i=1; i<=p; i++) printf(" Y%d", i); printf("\n"); printf("--------------------------------------------------------\n"); for (i=1; i3) printf(" %9.6f %9.6f %9.6f %9.6f %9.6f\n",x,t1[i],t2[i],t3[i],t4[i]); } printf("--------------------------------------------------------\n"); } /*************************************************************************** * SOLVING DIFFERENTIAL SYSTEMS WITH P VARIABLES OF ORDER 1 * * of type yi' = f(y1,y2,...,yn), i=1..n * * ------------------------------------------------------------------------ * * INPUTS: * * m number of points to display * * xi, xf begin, end values of variable x * * Yi table of begin values of functions at xi * * p number of independant variables * * fi finesse (number of intermediary points) * * * * OUTPUTS: * * t1,t2 real vectors storing the results for first two functions, * * y1 and y2. * * ------------------------------------------------------------------------ * * EXAMPLE: y1'=y2+y3-3y1, y2'=y1+y3-3y2, y3'=y1+y2-3y3 * * Exact solution : y1 = 1/3 (exp(-4x) + 2 exp(-x)) * * y2 = 1/3 (4exp(-4x) + 2 exp(-x)) * * y3 = 1/3 (-5exp(-4x)+ 2 exp(-x)) * ***************************************************************************/ void Eqdifp(double *t1,double *t2,double *t3,double *t4,double xi,double xf,double *Yi, int m,int p,int fi) { double h,x; double ta[10],tb[10],tc[10],td[10],y[10],z[10]; int i,j,k,ni; if (fi<1) return; h = (xf - xi) / fi / (m-1); p--; t1[1]=Yi[0]; t2[1]=Yi[1]; t3[1]=Yi[2]; t4[1]=Yi[3]; for (k=0; k