Back to index

wims  3.65+svn20090927
Romberg.java
Go to the documentation of this file.
00001 /*
00002  * Created on 05.11.2005
00003  *
00004  */
00005 package rene.zirkel.expression;
00006 
00007 import rene.zirkel.construction.ConstructionException;
00008 import rene.zirkel.objects.FunctionObject;
00009 
00010 public class Romberg
00011 {      
00014        private static double sumUp (FunctionObject F,
00015               double x, double h, int n)
00016        throws ConstructionException
00017        {   double sum=F.evaluateF(x);
00018          
00019            for (int i=1; i<=n; i++)
00020            {   sum=sum+F.evaluateF(x+i*h);
00021            }
00022            return sum;
00023        }
00024        
00034        static public double compute (FunctionObject F,
00035               double a, double b,
00036               int nstart, double eps, int maxiter)
00037        throws ConstructionException
00038        {                    
00039               // Ergebnisse der Trapezregel mit Schrittweite h/2^i
00040               double t[]=new double[maxiter];
00041               int n=nstart;
00042               
00043               // Anfangsschrittweite
00044               double h=(b-a)/n;
00045               // Berechne Trapezregel dazu und d[0]
00046               double tlast=t[0]=(F.evaluateF(a)+F.evaluateF(b)+
00047                                    2*sumUp(F,a+h,h,n-2))*h/2;
00048               double tlastbutone=tlast;
00049               // Bisheriges Ergebnis festhalten
00050               double old=t[0];
00051        
00052               // Halbiere Schrittweite bis Genauigkeit erreicht              
00053               for (int i=1; i<maxiter; i++)
00054               {      // Halbiere Schrittweite:
00055                      h=h/2;
00056                      n=n*2;
00057                      
00058                      // Berechne Trapezregel (unter Verwendung des
00059                      // letzten Ergebnisses:
00060                      t[i]=tlast/2+sumUp(F,a+h,2*h,n/2-1)*h;
00061                      tlastbutone=tlast;
00062                      tlast=t[i];
00063                      
00064                      // Update der t[i]:
00065                      double q=4;
00066                      for (int j=i-1; j>=0; j--)
00067                      {      t[j]=t[j+1]+(t[j+1]-t[j])/(q-1);
00068                             q=q*4;
00069                      }
00070                             
00071                      // Abbruch-Kriterium
00072                      double res=t[0];
00073                      if (Math.abs((res-old)/res)<eps) return res;
00074                      old=res;
00075               }
00076               
00077               // Bei ‹berschreiten der Interationsgrenze:
00078               return tlast;
00079        }
00080        
00081 }
00082 
00083