Check-ins Using Download Hex

Artifact 8b0feee3231260d0d1d6ebba036d462fe340dfaa621b5b4c1b766a209bd80324:

module numint; % compute a numerical integral

% Author: H. Melenk, 

     % Version March 1993

     % strategy for univariate integrals:
     %
     % 1. look for exact antiderivative, use that if it is 
     %          bounded in the interval,
     %
     % 2. try Chebyshev fit,
     %
     % 3. use M. Wulkow's adaptive multilevel quadrature.

% fixed a problem in finding the right accuracy 
% in the adaptive multilevel quadrature, W. Neun, August 1997

fluid '(!*noequiv accuracy!* singularities!*);
global '(iterations!* !*trnumeric);

symbolic procedure intrdeval u;
     % Interface function.
  begin scalar e,vars,y,p,singularities!*,imode,r,protfg!*,
          m,w;
    u := for each x in u collect reval x;
    u := accuracycontrol(u,3,20);
    e :=car u; u :=cdr u;
      % variables and interval's.
    % Allow for int(f,x,low,upp).
    if length u = 3 and (atom car u or not(caar u eq 'equal))
      then u := {{'equal,car u,'!*interval!* . cdr u}}
     else if eqcar(car u,'list) then
      u := for each x in cdar u collect reval x;
    for each x in u do
     <<if not eqcar(x,'equal) then typerr(x,"interval bounds");
       imode := 'infinity;
       if not !*rounded then setdmode('rounded,w:=!*rounded:=t);
       y:=revalnuminterval(caddr x,imode);
       if w then setdmode('rounded,w:=!*rounded:=nil);
       imode := t; % inf only for univariate cases.
       vars:=cadr x.vars;
       p:=(cadr x. car y. cadr y).p>>;
    m:=!*msg;
    if !*trnumeric then !*msg:=nil else protfg!* := t;  
    r:= intrd0(e,vars,p);
    erfg!* := nil; !*msg:=m;
    r:=reval r where dmode!* = '!:rd!:;
    return r;
  end;

put('num_int,'psopfn,'intrdeval);

symbolic procedure intrd0 (e,vars,p);
  if cdr vars then intrd2(e,vars,p) else
     % Univariate case.
   begin scalar fcn,x,lo,hi;
    fcn := e;
    x := car vars;
    lo := car cdar p; hi := cdr cdar p;
    return intrd1(fcn,x,lo,hi,0);
  end;

symbolic procedure intrd1(fcn,x,lo,hi,n);
 % Recursively decompose infinite intervals; map -inf to +inf;
 % transform integral over [a,inf] to integral over [0,1].
   begin scalar w,r1,r2,res;
    n := n+1;
    if !*trnumeric then
    <<writepri(n,'first); writepri(" level integrate :",nil);
      writepri(mkquote fcn,nil); writepri(" over ",nil);
      writepri(mkquote mkinterval(lo,hi),'last);
    >>;      
    if evalgreaterpx(lo,hi) then
    <<w:=lo;lo:=hi;hi:=w;fcn:={'minus,fcn}>>;

     % Map infinite intervals to finite ones.

     % case 1: [-inf,hi]

    if lo=minus!-infinity!* then
     if evalgreaterpx(hi,-1) then

     << % case 1.1: hi>-1 : [-inf,-1] + [-1,hi]
       r1 := intrd1(fcn,x,-1,hi,n);
       r2 := intrd1(subst({'minus,x},x,fcn),x,1,'infinity,n);
       res := reval {'plus,r1,r2};
     >> else

     << % case 1.2: hi <= -1
        res:= intrd1(subst({'minus,x},x,fcn),
                x,reval{'minus,lo},'infinity,n);
     >>;

     if res then goto fin;

     % case 2: [low .. inf] with -inf < low
    if hi='infinity then
     if evalgreaterpx(1,lo) then

     << % case 2.1: lo < 1: [lo,1] + [1,inf]
       r1 := intrd1(fcn,x,lo,1,n);
       r2 := intrd1(inttrans(fcn,1,x),x,0,1,n);
       res:= reval {'plus,r1,r2};
     >> else

     << % case 2.2: 1<=lo
        res := intrd1(inttrans(fcn,lo,x),x,0,1,n);
     >>;

     if res then goto fin;

     % case 3: finite interval
    res := intrd1a(fcn,x,lo,hi,{x.lo.hi});

 fin:
    if !*trnumeric then
    <<writepri(n,'first); writepri(" level integral :",nil);
      writepri(mkquote res,'last)>>;
    return res;
   end;         

symbolic procedure inttrans(u,a,x);
  reval inttrans1(u,a,x);

algebraic procedure inttrans1(u,a,x);
  % Transform function to integrate. The transformation
  % maps integral over [a,inf] to [0,1] when a > 0.
  % Here used only for a>1. 
  <<u:=sub(x=a/(1-!&t),u)*a /(a-!&t)^2;
    sub(!&t=x,u)>>;

symbolic procedure intrd1a(fcn,x,lo,hi,p);
   begin scalar u,w,acc,uu,loo,hii,oldmode,cbound;
     integer ord;
    if evalgreaterp(hi,lo) then <<loo:=lo; hii:=hi>>
      else <<loo:=hi; hii:=loo>>;
     % first attempt: look for bounded antiderivative.
    u := reval{'int,reval fcn,x} where dmode!*=nil,!*rounded=nil;
    oldmode:=switch!-mode!-rd nil;
    w:= errorset(
        {'boundseval,mkquote{u, {'equal,x,{'!*interval!*, loo, hii}}}},
                             nil,nil) where !*msg=nil;
    switch!-mode!-rd oldmode;
    
    if not smemq('int,u) and not errorp w then 
      <<
        if !*trnumeric then prin2t "Using bounded antiderivative";
        oldmode:=switch!-mode!-rd nil;
        u:=simp u;
        w:= prepsq subtrsq(subsq(u,{x . hi}),subsq(u,{x . lo}));
        switch!-mode!-rd oldmode;
        return w;
      >>;
    if !*trnumeric then prin2t "No bounded antiderivative found";

    % second attempt: look for a Chebyshev fit.
    w := nil;
    oldmode:=switch!-mode!-rd nil;
    acc := !:!:quotient(1,expt(10,accuracy!*));
    cbound := dm!: (acc /(hii - loo));
    ord := 20;
  chebloop:
    u := chebcoeffs(fcn,x,loo,hii,ord);
    uu := reverse u;
    if int!-chebconverges(uu,acc,cbound) then
    <<
      u:=chebint(u,nil,loo,hii); 
      w:=aeval{'difference,chebeval(u,nil,loo,hii,hi),
                           chebeval(u,nil,loo,hii,lo)};
    >>;
    if null w and ord<60 then <<ord:=ord+20; goto chebloop>>;
    switch!-mode!-rd oldmode;
    if w then
    <<if !*trnumeric then
       << prin2 "Using Chebyshev approximation of order "; 
          prin2t ord>>;
      return w;
    >>;
    if !*trnumeric then
    << 
       prin2t "No usable Chebyshev approximation found";
       prin2t "Starting adaptive multilevel quadrature";
    >>;
    return intrd2 (fcn,{x.lo.hi},p);
  end;

symbolic procedure int!-chebconverges(u,acc,ab);
  % Test the (reverse) Chebyshev coefficients u for convergence. 
  % acc: current accuracy (epsilon), ab: absolute upper bound for
  % trailing coefficient(s).
  begin scalar mx;
    mx := int!-chebmax(u,nil);
    dm!:(mx := mx*acc);
    return dm!:(abs car u < mx 
        and abs cadr u < mx and abs caddr u < mx
        and abs car u < ab and abs cadr u < ab 
        and abs caddr u < ab);
  end;

symbolic procedure  int!-chebmax(u,mx);
   if null u then mx else
   int!-chebmax(cdr u,if dm!:(w>mx) then w else mx) 
          where w=dm!: abs car u;
    

% ----------------- adaptive multilevel quadrature

symbolic procedure intrd2 (e,vars,p);
   begin scalar acc,r,oldmode,!*noequiv;
    vars := nil;
    oldmode:=switch!-mode!-rd nil;
    acc := !:!:quotient(1,expt(10,accuracy!*));
    e := reval e;
    r := if null cdr p then intrduni(e,p,acc) else
          intrdmulti(e,p,acc);
    switch!-mode!-rd oldmode;
    return r;
  end;

symbolic procedure intevaluate1(e,x,v);
  (if atom a then <<singularities!*:=v.singularities!*; 0>>
     else car a)where a=evaluate!*(e,list x,list v);

symbolic procedure intevaluate(e,x,v);
  (if atom a then <<singularities!*:=v.singularities!*; 0>>
     else car a)where a=evaluate!*(e,x,v);

symbolic procedure intrduni (e,p,acc);
   % univariate integral.
 dm!:
  begin scalar lo,hi,x,u,eps,i,il,int;
     integer n,k;
    x := car p; p:= cdr p;
   % evaluate the interval.
    lo := cadr x; hi := cddr x; x := car x;
    lo:=force!-to!-dm lo; hi:=force!-to!-dm hi;
    if hi=lo then return force!-to!-dm nil;
      % initial interval;
    il := list intrdmkinterval(e,x,
               (lo.intevaluate1(e,x,lo)),
               (hi.intevaluate1(e,x,hi)),1);
    int:=car lastpair car il;
  loop:
    k:=add1 k; n:=add1 n;
    if remainder(n,25)=0 then intrdmess(n,il);
   % divide the interval with the biggest local error;
    i:=car il; il:=cdr il; u:=intrdintersect(e,x,i);
    if null u then
        <<il:=i.il;
          intrdabortmsg(list car cadr i,list x,intcollecteps il);
           goto ready>>;
    for each q in u do il := intrdsortin(q,il);
    if k<5 then goto loop;

   % control variation of result
   %% if n=5 then   % W. Neun.
    int:=intcollectint il; % adjust accuracy
    k:=0; eps := intcollecteps il;
    if eps > (acc * abs int) then goto loop;
  ready:
    return intcollectint il;
  end;

symbolic procedure intrdabortmsg(pt,vars,eps);
   <<writepri(
"requested accuracy cannot be reached within iteration limit",'only);
     writepri("critical area of function near to ",'first);
     writepri(mkquote('list . for each u in pair(vars,pt)
       collect list('equal,car u,prepf cdr u)),'last);
     writepri("current error estimate: ",'first);
     writepri(mkquote prepf eps,'last);
   >>;

symbolic procedure intcollectint il;
 dm!: <<for each i in cdr il do r:= car lastpair i + r; r>>
          where r=car lastpair(car il);

symbolic procedure intcollecteps il;
 dm!: <<for each i in cdr il do r:= car i + r; r>>
          where r=car(car il);

symbolic procedure intrdsortin(u,l);
   % sort interval u into l such that l is sorted with descending car.
 dm!:( if null l then list u else
       if car u < caar l then car l . intrdsortin(u,cdr l) else
       u . l);

symbolic procedure intrdintersect(e,x,i);
  begin scalar d,plo,pmi,phi,lev;
    i:= cdr i;
    plo := car i; i := cdr i;
    pmi := car i; i := cdr i;
    phi := car i; i := cdr i;
    d   := car i; dm!:(d:=d/2);
    lev := cadr i +1;
    if lev>iterations!* then return nil;
    return list (intrdmkinterval(e,x,plo,pmi,lev) ,
                 intrdmkinterval(e,x,pmi,phi,lev) );
   end;

symbolic procedure intrdmkinterval(e,x,lo,hi,lev);
  dm!: begin scalar eps,it,is,mid,mi,flo,fhi,fmid,d,u;
    d := car hi- car lo;
    mid := (car hi+ car lo)/2;
    flo := cdr lo; fhi := cdr hi;
    fmid:=intevaluate1(e,x,mid);
    mi:=mid.fmid;
    if u:=intrdunisingular(lo,mi,hi) then
         <<flo:=car u; fmid:=cadr u; fhi:=caddr u>>;
    it := d*(flo+ 2fmid + fhi) / 4;         % order 1 integral;
    is := d*(4*flo + 16*fmid + 4*fhi)/24;   %simpson integral;
    eps:= abs(is-it);
                                            % interval record:
    return list (eps,                       % local error contrib
                 lo, mi, hi,                % the 3 points
                 d, lev, is);               % length and simpson.
  end;

symbolic procedure intrdunisingular(lo,mi,hi);
  % local extraploation for singular points.
  if singularities!* then
   begin scalar slo,smi,shi,fl,fm,fh,u;
    slo:=member(car lo,singularities!*);
    smi:=member(car mi,singularities!*);
    shi:=member(car hi,singularities!*);
    if not slo and not smi and not shi then return nil;
    if slo and smi and shi then rederr "too many singularities";
    fl:=cdr lo; fm:=cdr mi; fh:=cdr hi;
    dm!:( u := 2*(fl+fm+fh) );
    return list(if slo then u else fl,
                if smi then u else fm,
                if shi then u else fh);
  end;

symbolic procedure intrdmulti(e,p,acc);
  % multivariate integral.
 dm!:
  begin scalar vars,u,eps,i,il,int;
     integer n,k,dim;
    if smemq('infinity,p) then
      rederr "no support for infinity in multivariate integrals";
    dim:=length p;
    il:=intrdmkinitcube(e,p);
    vars := car il; il:= cdr il;
  loop:
    k:=add1 k; n:=add1 n;
    if remainder(n,25)=0 then intrdmess(n,il);
   % divide the simplex with the biggest local error;
    i:=car il; il:=cdr il; u:=intrdrefinecube(e,vars,i);
    if null u then
        <<il:=i.il;
          intrdabortmsg(car cadr i,vars,intcollecteps il);
           goto ready>>;
    for each q in u do il := intrdsortin(q,il);
    if k<5 then goto loop;

   % control variation of result
   % if n=5 then  % W. Neun.
    int:=intcollectint il; % adjust accuracy
    k:=0; eps := intcollecteps il;
    if eps> (acc * abs int) then goto loop;
  ready:
    return intcollectint il;
  end;

symbolic procedure intrdmkinitcube(e,p);
  % make initial element.
   begin scalar vol,vars,il,lo,hi,x,y;
    vol:=1;
    for each z in p do
    <<vars:=car z.vars;
      x:=force!-to!-dm cadr z;
      y:=force!-to!-dm cddr z;
      lo:=x.lo; hi:=y.hi;
      vol:= dm!:( abs(x-y)*vol );
    >>;
    lo:=lo.intevaluate(e,vars,lo);
    hi:=hi.intevaluate(e,vars,hi);
    il:=list intrdmkcube(e,vars,lo,hi,vol,1);
    return vars.il;
  end;

symbolic procedure intrdrefinecube(e,vars,i);
 dm!:  % divide cube into 2.
  begin scalar plo,phi,lo,hi,nlo,nhi,vol,x,y,xnew;
     integer m,lev;
    i:=cdr i;
    plo:=car i; i:=cdr i; phi:=car i; i:=cdr i;
    vol:= car i / 2; lev := add1 cadr i;
    if lev > iterations!* then return nil;
    lo:=car plo; hi:=car phi;
      % select longest edge of interval;
    m := 1; x := car hi-car lo;
    for j:=2:length lo do
      if x<(y:=(nth(hi,j)-nth(lo,j))) then<<m:=j;x:=y>>;
    nlo := append(lo,nil); nhi := append(hi,nil);
    xnew:=(nth(hi,m) + nth(lo,m)) /2 ;
    nth(nlo,m):=xnew; nth(nhi,m):=xnew;
    nlo := nlo.intevaluate(e,vars,nlo);
    nhi := nhi.intevaluate(e,vars,nhi);
    return list(
      intrdmkcube(e,vars,plo,nhi,vol,lev),
      intrdmkcube(e,vars,nlo,phi,vol,lev));
  end;

symbolic procedure intrdmkcube(e,vars,lo,hi,vol,lev);
   % make a fresh cube record.
  dm!: begin scalar u,eps,it,is,mid,fmi,flo,fhi;
    flo:= cdr lo; fhi:= cdr hi;
    mid:=list!+list(car lo,car hi); mid:=scal!*list(1/2,mid);
    fmi:= intevaluate(e,vars,mid);
    if u:=intrdunisingular(lo,mid.fmi,hi) then
         <<flo:=car u; fmi:=cadr u; fhi:=caddr u>>;
    is:=flo+fhi;
    it:=is*vol/2; is:=(2*fmi+is)*vol/4;
    eps := abs(is-it);
    return list(eps,lo,hi,vol,lev,is);
  end;

symbolic procedure intrdmess(n,il);
  if !*trnumeric then
   <<writepri(2*n,'first); writepri(" intervals, integral=",nil);
     writepri(mkquote prepf intcollectint il,nil);
     writepri(", error estimate=",nil);
     writepri(mkquote prepf intcollecteps il,nil);
   >>;

symbolic procedure prinxt l;
  begin integer i;
    for each x in l do
     if not eqcar(x,'!:rd!:) then prin2 x else
     if atom cdr x then prin2 x else
       <<prin2 (float cadr x *expt(10.0,cddr x));
         tab (i:=20+i)>>;
     terpri();
   end;

endmodule;

end;
REDUCE Historical
[ REDUCE Sourceforge Project | Historical SVN Repository | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]