/************************************************************************* * Solve an ordinary system of first order differential equations (N<=10) * * with initial conditions using a Runge-Kutta integration method. * * ---------------------------------------------------------------------- * * SAMPLE RUN: * * Integrate first order DE system: * * y1' = y2 * * y2' = y3 * * y3' = y4 * * y4' = y5 * * y5' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3) * * from x=0 to x=2 with initial conditions: y(0)..y(4) = 1 * * (25 integration steps, finesse=20). * * * * X Y1 <------------------------> YN estimated * * ---------------------------------------------------------------------- * * 0.0000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 * * 0.0800 1.0832871 1.0832862 1.0832443 1.0816193 1.0382909 * * 0.1600 1.1735104 1.1734962 1.1731252 1.1657018 1.0606156 * * 0.2400 1.2712454 1.2711671 1.2697821 1.2507235 1.0606138 * * 0.3200 1.3771106 1.3768421 1.3732111 1.3346056 1.0307873 * * 0.4000 1.4917679 1.4910565 1.4832168 1.4146287 0.9626611 * * 0.4800 1.6159210 1.6143206 1.5993581 1.4873681 0.8471623 * * 0.5600 1.7503129 1.7470979 1.7208913 1.5486703 0.6752862 * * 0.6400 1.8957208 1.8897779 1.8467119 1.5936970 0.4391094 * * 0.7200 2.0529492 2.0426461 1.9753055 1.6170637 0.1331626 * * 0.8000 2.2228199 2.2058489 2.1047132 1.6131014 -0.2438996 * * 0.8800 2.4061600 2.3793580 2.2325243 1.5762584 -0.6875083 * * 0.9600 2.6037859 2.5629350 2.3559061 1.5016398 -1.1855881 * * 1.0400 2.8164856 2.7561004 2.4716817 1.3856519 -1.7176006 * * 1.1200 3.0449981 2.9581109 2.5764618 1.2266816 -2.2547153 * * 1.2000 3.2899923 3.1679497 2.6668279 1.0257030 -2.7615081 * * 1.2800 3.5520448 3.3843326 2.7395571 0.7866810 -3.1993222 * * 1.3600 3.8316195 3.6057337 2.7918674 0.5166480 -3.5309960 * * 1.4400 4.1290502 3.8304300 2.8216525 0.2253716 -3.7261918 * * 1.5200 4.4445267 4.0565635 2.8276730 -0.0753912 -3.7662105 * * 1.6000 4.7780875 4.2822163 2.8096752 -0.3729657 -3.6471428 * * 1.6800 5.1296178 4.5051904 2.7684146 -0.6549877 -3.3805446 * * 1.7600 5.4988544 4.7245867 2.7055842 -0.9105780 -2.9914571 * * 1.8400 5.8853965 4.9378740 2.6236555 -1.1312684 -2.5142924 * * 1.9200 6.2887215 5.1439429 2.5256617 -1.3115454 -1.9876112 * * 2.0000 6.7082039 5.3416408 2.4149534 -1.4489720 -1.4489720 * * * * ---------------------------------------------------------------------- * * 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. * * (www.jpmoreau.fr) * *************************************************************************/ #include #include #define SIZE 10 //maximum number of DEs double X,Y[SIZE],Y1[SIZE]; double xl,x0,h; int finesse,i,k,kl,l,n; // DE system to integrate (n=5) void F(double x, double *y, double *yp) { yp[0] = y[1]; yp[1] = y[2]; yp[2] = y[3]; yp[3] = y[4]; yp[4] = (45.0 * y[2] * y[3] * y[4] - 40.0 * y[3] * y[3] * y[3]) / (9.0 * y[2] * y[2]); } //Integrate first order DE sytem from x to x+h using Runge-Kutta void RK4n(int n, double x, double *y, double h, double *y1) { double c1[SIZE],c2[SIZE],c3[SIZE],c4[SIZE],yy[SIZE],h2; int i; F(x,y,c1); h2=h/2.0; for (i=0; i YN estimated\n"); printf(" ----------------------------------------------------------------------\n"); // integration loop RK4n(n,X,Y,h,Y1); printf("%9.4f", X); for (i=0; i