I have been thinking about this & I don't see a better solution than
Raul's. @Raul: think about putting something in NuVoc explaining this.
I thought at first: should the anonymous verb created by
verb adv_locale_
automatically be executed in (locale)? That would solve the immediate
problem.
But it leaves us with the responsibility of defining a locale for every
anonymous verb. What locale should we assign to:
(V0 V1_locale_ V2)
(V0_locale0_ V1_locale1_ V2_locale2_)
(V0_locale0_ V1)
?
I haven't been able to come up with something I would be willing to
suggest to Roger & Ken.
Unless someone does, I guess we will leave it as is. The locale suffix
gives an easy way to specify the locale of a verb - 98% of the cases -
and for the rest you need Raul's Device.
Henry Rich
On 4/4/2017 4:03 AM, Raul Miller wrote:
I am not sure how to test this (I do not know what values to use for f
a and b in 'f qromb a b' and do not feel like studying your code
enough to deduce such values), but I think this is what you are asking
for:
---------- 8< ------ cut here ------ 8< ----------
coclass'nr'
qromb=:1 :0
qromb_implementation (coname '')
)
qromb_implementation=:2 :0
NB.Returns the integral of the function from a to b. NR P.166 Ch.4.3
NB.Uses Romberg's method of order 2*K, where, e.g. K=2 is Simpson's rule.
NB.u function
NB.y a,b,eps integration bounds(a,b), and accuracy(eps,
default 1e_10)
NB.eps is the fraction error from extrapolation error estimate
NB.PolyInterpX stores successive trapezoidal relative stepsizes
NB.PolyInterpY stores their approximations
NB.K is the number of points used in the extrapolation
cocurrent v
ab=.2{.y [ eps=.2{y,1e_10
PolyInterpY=:20#0 [ PolyInterpX=:21#0 [ PolyInterpM=:K=.5
PolyInterpX=:PolyInterpX 0}~1
TrapzdNextN=:0
j=.1 while.j<:#PolyInterpX do.
PolyInterpY=:PolyInterpY(j-1)}~u trapzdNext ab
if.j>:K do.'ss dy'=.0 polyInterpRawinterp j-K
if.(|dy)<:eps*|ss do.ss return.end.end.
NB.This is key. The factor 0.25 allows h^2 extrapolation. See NR
equation 4.2.1.
PolyInterpX=:PolyInterpX j}~0.25*PolyInterpX{~j-1
j=.>:j end.
'Too many steps in routine qromb'assert 0
)
polyInterpRawinterp=:4 :0
NB.Polynomial interpolation. NR P.119 Ch.3.2
NB.x the point of interpolation
NB.y j subrange j+i.PolyInterpM is used for the interpolation
NB.Must initialize
NB.PolyInterpM=:5
NB.PolyInterpX=:21#0
NB.PolyInterpY=:20#0
j=.y
dif=.|x-j{PolyInterpX
i=.0 while.i<PolyInterpM do.
if.dif>dift=.|x-PolyInterpX{~j+i do.ns=.i [ dif=.dift end.
i=.>:i end.
d=.c=.PolyInterpY ];.0~ j,:PolyInterpM
ns=.<:ns [ y=.PolyInterpY{~j+ns
m=.1 while.m<PolyInterpM do.
i=.0 while.i<PolyInterpM-m do.
ho=.x-~PolyInterpX{~j+i
hp=.x-~PolyInterpX{~j+i+m
w=.(c{~i+1)-i{d
'PolyInterp error'assert 0~:den=.ho-hp
den=.w%den
d=.d i}~hp*den
c=.c i}~ho*den
i=.>:i end.
if.(PolyInterpM-m)>2*ns+1 do.dy=.c{~ns+1
else.ns=.<:ns [ dy=.ns{d
end.
y=.y+dy
m=.>:m end.
y,dy
)
trapzdNext=:1 :0
trapzdNext_implementation (coname '')
)
trapzdNext_implementation=:2 :0
NB.Returns the nth stage of refinement of the extended trapezoidal
rule. NR P.163 Ch.4.2
NB.u function must accept list
NB.y a,b range
NB.Must initialize TrapzdNextN=:0 before using.
cocurrent v
ba=.-~/y
TrapzdNextN=:>:TrapzdNextN
if.1=TrapzdNextN do.TrapzdNextS=:-:ba*+/u y return.
else.TrapzdNextS=:-:TrapzdNextS+ba*t%~+/u({.y)+(0.5+i.t)*ba%t=.2^TrapzdNextN-2
return.
end.
)
---------- 8< ------ cut here ------ 8< ----------
This is a rather bulky example, if there are problems with this
approach it might be better to define a more concise (and complete -
with a test case which illustrates the problem) example?
Thanks,
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
