/* Versions 1. Romberg.ox, November 2000, Anders Karlstrom 2. Romberg_.ox, April 19. 2003 Anders Karlsirom and Edward Morey This program creates a number of functions that are used by the program ExpectedCV.ox. Specifically, the functions created in this program is used by that program This program cannot be run by itself (it uses a number of functions defined in ExpectedCV.ox), so should be viewed as part of that program. It defines the functions trapzd(), polint(), qromb() which is used to solve the integral of the function "integrand" in ExpectedCV.ox. The method is Romberg integration as outlined in Numerical Recipes. The function to be integrated should return a vector. The user does not have to modify this file. */ trapzd(const func,const a,const b,const s_ut, const n) // "func" is the function to be integrated { decl m_s; decl it,j; decl del,sum,tnm,x; m_s= s_static; if (n==1) m_s = 0.5*(b-a) .* (func(a)+func(b)); else { it=2^(n-2); tnm=it; del=(b-a)/tnm; x=a+0.5*del; sum=0; for (j=1;j<=it;j++) { sum=sum+func(x); x=x+del; } m_s=0.5*(m_s+(b-a) .* sum ./ tnm); } s_ut[0] = m_s'; s_static=m_s; return 0; } polint(const xa, const ya, const n, const x, const y_ut, const dy_ut) { // in: xa and ya, degree n, approximation at x, out: y_ut, dy_ut decl NMAX; decl y=y_ut,dy=dy_ut; NMAX=10; decl i,m,ns; decl den,dif,dift,ho,hp,w,c=zeros(NMAX,1),d=zeros(NMAX,1) ; ns=1; dif=abs(x-xa[1][]); for (i=1;i<=n;i++) { dift=abs(x-xa[i][]); if (diftK) { h_tmp = h[j-KM-1:]; s_tmp=s[j-KM-1:][]; polint(h_tmp*ones(1,mNOBS),s_tmp,K,zeros(1,mNOBS),&m_ss,&dss); diff=m_ss-mLastResult; maxerr = limits(dss'); // print("Correction in this iteration:", limits((diff)')); // println("Elapsed time:",timespan(tid)); // println("--------------------------------------"); mLastResult=m_ss; if (abs(diff)