RosettaCodeData/Task/Numerical-integration/Phix/numerical-integration.phix

75 lines
2.3 KiB
Plaintext

function rect_left(integer rid, atom x, atom h)
if atom(h) then end if -- suppress warning
return call_func(rid,{x})
end function
function rect_mid(integer rid, atom x, atom h)
return call_func(rid,{x+h/2})
end function
function rect_right(integer rid, atom x, atom h)
return call_func(rid,{x+h})
end function
function trapezium(integer rid, atom x, atom h)
return (call_func(rid,{x})+call_func(rid,{x+h}))/2
end function
function simpson(integer rid, atom x, atom h)
return (call_func(rid,{x})+4*call_func(rid,{x+h/2})+call_func(rid,{x+h}))/6
end function
function cubed(atom x)
return power(x,3)
end function
function recip(atom x)
return 1/x
end function
function ident(atom x)
return x
end function
function integrate(integer m_id, integer f_id, atom a, atom b, integer steps)
atom accum = 0,
h = (b-a)/steps
for i=0 to steps-1 do
accum += call_func(m_id,{f_id,a+h*i,h})
end for
return h*accum
end function
function smartp(atom N)
string res
if N=floor(N) then return sprintf("%d",N) end if
res = sprintf("%12f",round(N,1000000))
if find('.',res) then
res = trim_tail(res,"0")
res = trim_tail(res,".")
end if
return res
end function
procedure test(sequence tests)
string name
atom a, b, steps, rid
printf(1,"Function Range Iterations L-Rect M-Rect R-Rect Trapeze Simpson\n")
for i=1 to length(tests) do
{name,a,b,steps,rid} = tests[i]
printf(1," %-5s %6d - %-5d %10d %12s %12s %12s %12s %12s\n",{name,a,b,steps,
smartp(integrate(routine_id("rect_left"), rid,a,b,steps)),
smartp(integrate(routine_id("rect_mid"), rid,a,b,steps)),
smartp(integrate(routine_id("rect_right"), rid,a,b,steps)),
smartp(integrate(routine_id("trapezium"), rid,a,b,steps)),
smartp(integrate(routine_id("simpson"), rid,a,b,steps))})
end for
end procedure
constant tests = {{"x^3", 0, 1, 100, routine_id("cubed")},
{"1/x", 1, 100, 1000, routine_id("recip")},
{"x", 0, 5000, 5000000, routine_id("ident")},
{"x", 0, 6000, 6000000, routine_id("ident")}}
test(tests)