Rosetta Code

Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

Page Discussion

Numerical integration/Gauss-Legendre Quadrature (view source)

Revision as of 18:24, 14 March 2022

15,548 bytes added ,  2 years ago
m
→‎{{header|Phix}}: syntax coloured
(→‎{{header|Wren}}: Now uses new core library method.)
m (→‎{{header|Phix}}: syntax coloured)
Line 1,777:
=={{header|Phix}}==
{{trans|Lua}}
<!--<lang Phix>integer order = 0(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">order</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
sequence legendreRoots = {},
legendreWeights = {}
<span style="color: #004080;">sequence</span> <span style="color: #000000;">legendreRoots</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{},</span>
<span style="color: #000000;">legendreWeights</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
function legendre(integer term, atom z)
if term=0 then
<span style="color: #008080;">function</span> <span style="color: #000000;">legendre</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
return 1
<span style="color: #008080;">if</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
elsif term=1 then
<span style="color: #008080;">return</span> <span style="color: #000000;">1</span>
return z
<span style="color: #008080;">elsif</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span>
else
<span style="color: #008080;">return</span> <span style="color: #000000;">z</span>
return ((2*term-1)*z*legendre(term-1,z)-(term-1)*legendre(term-2,z))/term
<span style="color: #008080;">else</span>
end if
<span style="color: #008080;">return</span> <span style="color: #0000FF;">((</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">term</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">z</span><span style="color: #0000FF;">*</span><span style="color: #000000;">legendre</span><span style="color: #0000FF;">(</span><span style="color: #000000;">term</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)-(</span><span style="color: #000000;">term</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">legendre</span><span style="color: #0000FF;">(</span><span style="color: #000000;">term</span><span style="color: #0000FF;">-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">))/</span><span style="color: #000000;">term</span>
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
function legendreDerivative(integer term, atom z)
if term=0
<span style="color: #008080;">function</span> <span style="color: #000000;">legendreDerivative</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
or term=1 then
<span style="color: #008080;">if</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span>
return term
<span style="color: #008080;">or</span> <span style="color: #000000;">term</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span>
end if
<span style="color: #008080;">return</span> <span style="color: #000000;">term</span>
return (term*(z*legendre(term,z)-legendre(term-1,z)))/(z*z-1)
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
end function
<span style="color: #008080;">return</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">term</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">*</span><span style="color: #000000;">legendre</span><span style="color: #0000FF;">(</span><span style="color: #000000;">term</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)-</span><span style="color: #000000;">legendre</span><span style="color: #0000FF;">(</span><span style="color: #000000;">term</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)))/(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">*</span><span style="color: #000000;">z</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
procedure getLegendreRoots()
legendreRoots = {}
<span style="color: #008080;">procedure</span> <span style="color: #000000;">getLegendreRoots</span><span style="color: #0000FF;">()</span>
for index=1 to order do
<span style="color: #000000;">legendreRoots</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
atom y = cos(PI*(index-0.25)/(order+0.5))
<span style="color: #008080;">for</span> <span style="color: #000000;">index</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">order</span> <span style="color: #008080;">do</span>
while 1 do
<span style="color: #004080;">atom</span> <span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">cos</span><span style="color: #0000FF;">(</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">index</span><span style="color: #0000FF;">-</span><span style="color: #000000;">0.25</span><span style="color: #0000FF;">)/(</span><span style="color: #000000;">order</span><span style="color: #0000FF;">+</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">))</span>
atom y1 = y
<span style="color: #008080;">while</span> <span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
y -= legendre(order,y)/legendreDerivative(order,y)
<span style="color: #004080;">atom</span> <span style="color: #000000;">y1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">y</span>
if abs(y-y1)<2e-16 then exit end if
<span style="color: #000000;">y</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">legendre</span><span style="color: #0000FF;">(</span><span style="color: #000000;">order</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">legendreDerivative</span><span style="color: #0000FF;">(</span><span style="color: #000000;">order</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">)</span>
end while
<span style="color: #008080;">if</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">y</span><span style="color: #0000FF;">-</span><span style="color: #000000;">y1</span><span style="color: #0000FF;">)<</span><span style="color: #000000;">2e-16</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
legendreRoots &= y
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
end for
<span style="color: #000000;">legendreRoots</span> <span style="color: #0000FF;">&=</span> <span style="color: #000000;">y</span>
end procedure
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
procedure getLegendreWeights()
legendreWeights = {}
<span style="color: #008080;">procedure</span> <span style="color: #000000;">getLegendreWeights</span><span style="color: #0000FF;">()</span>
for index=1 to order do
<span style="color: #000000;">legendreWeights</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
atom lri = legendreRoots[index],
<span style="color: #008080;">for</span> <span style="color: #000000;">index</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">order</span> <span style="color: #008080;">do</span>
diff = legendreDerivative(order,lri),
<span style="color: #004080;">atom</span> <span style="color: #000000;">lri</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">legendreRoots</span><span style="color: #0000FF;">[</span><span style="color: #000000;">index</span><span style="color: #0000FF;">],</span>
weight = 2 / ((1-power(lri,2))*power(diff,2))
<span style="color: #000000;">diff</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">legendreDerivative</span><span style="color: #0000FF;">(</span><span style="color: #000000;">order</span><span style="color: #0000FF;">,</span><span style="color: #000000;">lri</span><span style="color: #0000FF;">),</span>
legendreWeights &= weight
<span style="color: #000000;">weight</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">2</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">((</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lri</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">))*</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">diff</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">))</span>
end for
<span style="color: #000000;">legendreWeights</span> <span style="color: #0000FF;">&=</span> <span style="color: #000000;">weight</span>
end procedure
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
function gaussLegendreQuadrature(integer f, lowerLimit, upperLimit, n)
order = n
<span style="color: #008080;">function</span> <span style="color: #000000;">gaussLegendreQuadrature</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lowerLimit</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">upperLimit</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">order</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">n</span>
getLegendreRoots()
getLegendreWeights()
<span style="color: #000000;">getLegendreRoots</span><span style="color: #0000FF;">()</span>
<span style="color: #000000;">getLegendreWeights</span><span style="color: #0000FF;">()</span>
atom c1 = (upperLimit - lowerLimit) / 2
atom c2 = (upperLimit + lowerLimit) / 2
<span style="color: #004080;">atom</span> <span style="color: #000000;">c1</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">upperLimit</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">lowerLimit</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">2</span>
atom s = 0
<span style="color: #004080;">atom</span> <span style="color: #000000;">c2</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">upperLimit</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">lowerLimit</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">2</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
for i = 1 to order do
s += legendreWeights[i] * call_func(f,{c1 * legendreRoots[i] + c2})
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">order</span> <span style="color: #008080;">do</span>
end for
<span style="color: #000000;">s</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">legendreWeights</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">c1</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">legendreRoots</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">c2</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
return c1 * s
end function
<span style="color: #008080;">return</span> <span style="color: #000000;">c1</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">s</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
include pmaths.e -- exp()
constant r_exp = routine_id("exp")
<span style="color: #004080;">string</span> <span style="color: #000000;">fmt</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">machine_bits</span><span style="color: #0000FF;">()=</span><span style="color: #000000;">32</span><span style="color: #0000FF;">?</span><span style="color: #008000;">"%.13f"</span><span style="color: #0000FF;">:</span><span style="color: #008000;">"%.14f"</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">5</span> <span style="color: #008080;">to</span> <span style="color: #000000;">11</span> <span style="color: #008080;">by</span> <span style="color: #000000;">6</span> <span style="color: #008080;">do</span>
string fmt = iff(machine_bits()=32?"%.13f":"%.14f")
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">gaussLegendreQuadrature</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">exp</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">)})</span>
string res
<span style="color: #008080;">if</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">5</span> <span style="color: #008080;">then</span>
for i=5 to 11 by 6 do
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"roots:"</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">legendreRoots</span>
res = sprintf(fmt,{gaussLegendreQuadrature(r_exp, -3, 3, i)})
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"weights:"</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">legendreWeights</span>
if i=5 then
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
puts(1,"roots:") ?legendreRoots
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Gauss-Legendre %2d-point quadrature for exp over [-3..3] = %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">order</span><span style="color: #0000FF;">,</span><span style="color: #000000;">res</span><span style="color: #0000FF;">})</span>
puts(1,"weights:") ?legendreWeights
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end if
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">3</span><span style="color: #0000FF;">)-</span><span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">)})</span>
printf(1,"Gauss-Legendre %2d-point quadrature for exp over [-3..3] = %s\n",{order,res})
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" compared to actual = %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">res</span><span style="color: #0000FF;">})</span>
end for
<!--</lang>-->
res = sprintf(fmt,{exp(3)-exp(-3)})
printf(1," compared to actual = %s\n",{res})</lang>
{{out}}
<pre>
Petelomax
7,795

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/183831"