Numerical integration: Difference between revisions

Line 3,955:

=={{header|Phix}}==

<lang Phix>function rect_left(integer rid, atom x, atom h)

function rect_left(integer rid, atom x, atom /*h*/)

if atom(h) then end if -- suppress warning

return rid(x)

return call_func(rid,{x})

end function

end function

function rect_mid(integer rid, atom x, atom h)

function rect_mid(integer rid, atom x, atom h)

return rid(x+h/2)

return call_func(rid,{x+h/2})

end function

end function

function rect_right(integer rid, atom x, atom h)

function rect_right(integer rid, atom x, atom h)

return rid(x+h)

return call_func(rid,{x+h})

end function

end function

function trapezium(integer rid, atom x, atom h)

function trapezium(integer rid, atom x, atom h)

return (rid(x)+rid(x+h))/2

return (call_func(rid,{x})+call_func(rid,{x+h}))/2

end function

end function

function simpson(integer rid, atom x, atom h)

function simpson(integer rid, atom x, atom h)

return (rid(x)+4*rid(x+h/2)+rid(x+h))/6

return (call_func(rid,{x})+4*call_func(rid,{x+h/2})+call_func(rid,{x+h}))/6

end function

end function

function cubed(atom x)

function cubed(atom x)

return power(x,3)

return power(x,3)

end function

end function

function recip(atom x)

function recip(atom x)

return 1/x

return 1/x

end function

end function

function ident(atom x)

function ident(atom x)

return x

return x

end function

end function

function integrate(integer m_id, integer f_id, atom a, atom b, integer steps)

function integrate(integer m_id, integer f_id, atom a, atom b, integer steps)

atom accum = 0,

atom accum = 0,

h = (b-a)/steps

h = (b-a)/steps

for i=0 to steps-1 do

for i=0 to steps-1 do

accum += m_id(f_id,a+h*i,h)

accum += call_func(m_id,{f_id,a+h*i,h})

end for

end for

return h*accum

return h*accum

end function

end function

function smartp(atom N)

function smartp(atom N)

if N=floor(N) then return sprintf("%d",N) end if

string res

string res = sprintf("%12f",round(N,1000000))

if N=floor(N) then return sprintf("%d",N) end if

if find('.',res) then

res = sprintf("%12f",round(N,1000000))

res = trim_tail(res,"0")

if find('.',res) then

res = trim_tail(res,".")

res = trim_tail(res,"0")

end if

res = trim_tail(res,".")

return res

end if

end function

return res

end function

procedure test(sequence tests)

string name

procedure test(sequence tests)

atom a, b, steps, rid

string name

printf(1,"Function Range Iterations L-Rect M-Rect R-Rect Trapeze Simpson\n")

atom a, b, steps, rid

for i=1 to length(tests) do

printf(1,"Function Range Iterations L-Rect M-Rect R-Rect Trapeze Simpson\n")

{name,a,b,steps,rid} = tests[i]

for i=1 to length(tests) do

printf(1," %-5s %6d - %-5d %10d %12s %12s %12s %12s %12s\n",{name,a,b,steps,

{name,a,b,steps,rid} = tests[i]

smartp(integrate(rect_left, rid,a,b,steps)),

printf(1," %-5s %6d - %-5d %10d %12s %12s %12s %12s %12s\n",{name,a,b,steps,

smartp(integrate(rect_mid, rid,a,b,steps)),

smartp(integrate(routine_id("rect_left"), rid,a,b,steps)),

smartp(integrate(rect_right,rid,a,b,steps)),

smartp(integrate(routine_id("rect_mid"), rid,a,b,steps)),

smartp(integrate(trapezium, rid,a,b,steps)),

smartp(integrate(routine_id("rect_right"), rid,a,b,steps)),

smartp(integrate(simpson, rid,a,b,steps))})

smartp(integrate(routine_id("trapezium"), rid,a,b,steps)),

end for

smartp(integrate(routine_id("simpson"), rid,a,b,steps))})

end procedure

end for

end procedure

constant tests = {{"x^3", 0, 1, 100, cubed},

{"1/x", 1, 100, 1000, recip},

constant tests = {{"x^3", 0, 1, 100, routine_id("cubed")},

{"x", 0, 5000, 5000000, ident},

{"1/x", 1, 100, 1000, routine_id("recip")},

{"x", 0, 6000, 6000000, ident}}

{"x", 0, 5000, 5000000, routine_id("ident")},

{"x", 0, 6000, 6000000, routine_id("ident")}}

test(tests)

~~test(tests)~~</lang>

<pre>