The "C + +" Implementation of the Dragon backlog algorithm

1. The formula of the re-trapezoid method and the recursive

The method of complex trapezoid is an effective way to improve the precision of quadrature formula. Divide the [a, b] interval n, step h = (b-a)/n, point xk = a + kh. The formula of complex quadrature is to calculate the normal trapezoidal method between each cell in this n equal section, and then sum up the small interval of n. The formula is as follows:

When using the complex trapezoidal method integral, this process can be recursive, to more convenient use of computer implementation. Set the integral interval [a, b], the interval n equal, then the sub-point total of n+1, using the re-trapezoidal integral to obtain TN. To make two points, the results of the dichotomy are recorded as t2n, there are:

2. Romberg Integral formula

Romberg integral is actually an algorithm to improve the convergence rate. Due to the reduction of the error after the step length of the re-trapezoid method is reduced to:

Organized by:

According to this idea, the convergent slow trapezoid sequence tn is processed into a fast convergent Romberg value Sequence RN, this is the Romberg algorithm, the processing algorithm flow is as follows:

Realize:

1#include <stdio.h>2#include <math.h>3#include <iostream>4#include <cstdio>5 using namespacestd;6 intrk=0;7 inttk=0;8 DoubleFxDoubleX//integrable function9 {Ten     //if (x==0.0) return 1.0; One     return 3*x*x*x+2*x*x+1+sin (x); A } - DoubleGetreal (DoubleADoubleb) { -     DoubleR1 =3.0/4.0* B*b*b*b +2.0/3.0*b*b*b + B-cos (b); the     Doubler2 =3.0/4.0* A*a*a*a +2.0/3.0*a*a*a + A-cos (a); -     returnR1-R2; - } - DoubleGetS (DoubleADoubleBDoubleh) + { -     Doubleres=0.0; +      for(Doublei=a+h/2.0; i<b; i+=h) { Ares+=FX (i);  at     } -          -     returnRes; - } - DoubleRomberg (DoubleADoubleBDoublee) - { in     intk=1; -     Doublet1,t2,s1,s2,c1,c2,r1,r2; to     Doubleh=b-A; +     Doubles; -T1= (FX (a) +fx (b)) *h/2.0; the     intCounter=0; *      while(1) $     {Panax Notoginsengrk++; -counter++; thes=GetS (a,b,h); +T2= (t1+h*s)/2.0; AS2= (4.0*T2-T1)/3.0; theH/=2.0; +t1=T2; -s1=S2; $c1=C2; $r1=R2; -         if(k==1) -         { thek++; -             Continue;Wuyi         } theC2= (16.0*S2-S1)/15.0; -         if(k==2) Wu         { -k++; About             Continue; $         } -R2= (64.0*C2-C1)/63.0; -         if(k==3) -         { Ak++; +             Continue; the         } -         if(Fabs (R1-R2) <e| | counter>= -) Break; $     } the     returnR2; the } the DoubleTn (DoubleADoubleBDoublee) the { -     Doublet1,t2; in     Doubleh=b-A; theT1= (FX (a) +fx (b)) *h/2.0; the      while(1) About     { thetk++; the         Doubles=GetS (a,b,h); theT2= (t1+h*s)/2.0; +         if(Fabs (T2-T1) <e) Break; -H/=2.0; thet1=T2;Bayi     } the     returnT2; the } - intMain () - { the     Doublea,b,e; theprintf"input integral limit and accuracy: a b E:"); the     //input interval [a, b], and precision e thescanf"%LF%LF%LF",&a,&b,&e); -     Doublet=Romberg (a,b,e); the     //the results of the Romberg algorithm and trapezoid method are respectively output and the corresponding two-point number theprintf"\nromberg: Integral value:%.7LF--Two minutes:%d\n", t,rk); thet=Tn (a,b,e);94printf"Tn: Integral value:%.7LF--Two minutes:%d\n", T,TK); the     Doubletf =Getreal (A, b); theprintf"REAL:%.7LF", TF); the     return 0;98}

The "C + +" Implementation of the Dragon backlog algorithm