Numerical integration/Gauss-Legendre Quadrature: Difference between revisions
Content added Content deleted
Walterpachl (talk | contribs) (→{{header|PL/I}}: correct one declaration for the better!?!) |
Walterpachl (talk | contribs) (REXX added (translated from PL/I)) |
||
Line 989: | Line 989: | ||
</lang> |
</lang> |
||
[[File:gauss.png]] |
[[File:gauss.png]] |
||
=={{header|REXX}}== |
|||
<lang rexx>/*--------------------------------------------------------------------- |
|||
* 31.10.2013 Walter Pachl Translation from PL/I |
|||
*--------------------------------------------------------------------*/ |
|||
prec=60 |
|||
Numeric Digits prec |
|||
epsilon=1/10**prec |
|||
pi=3.141592653589793238462643383279502884197169399375105820974944592307 |
|||
exact = exp(3,prec)-exp(-3,prec) |
|||
Do n = 1 To 20 |
|||
a = -3; b = 3 |
|||
r.=0 |
|||
call gaussquad |
|||
sum=0 |
|||
Do j=1 To n |
|||
sum=sum + r.2.j * exp((a+b)/2+r.1.j*(b-a)/2,prec) |
|||
End |
|||
z = (b-a)/2 * sum |
|||
Say right(n,2) format(z,2,40) format(z-exact,2,4,,0) |
|||
End |
|||
Say ' ' exact '(exact)' |
|||
Exit |
|||
gaussquad: |
|||
p0.0=1; p0.1=1 |
|||
p1.0=2; p1.1=1; p1.2=0 |
|||
Do k = 2 To n |
|||
tmp.0=p1.0+1 |
|||
Do L = 1 To p1.0 |
|||
tmp.l = p1.l |
|||
End |
|||
tmp.l=0 |
|||
tmp2.0=p0.0+2 |
|||
tmp2.1=0 |
|||
tmp2.2=0 |
|||
Do L = 1 To p0.0 |
|||
l2=l+2 |
|||
tmp2.l2=p0.l |
|||
End |
|||
Do j=1 To tmp.0 |
|||
tmp.j = ((2*k-1)*tmp.j - (k-1)*tmp2.j)/k |
|||
End |
|||
p0.0=p1.0 |
|||
Do j=1 To p0.0 |
|||
p0.j = p1.j |
|||
End |
|||
p1.0=tmp.0 |
|||
Do j=1 To p1.0 |
|||
p1.j=tmp.j |
|||
End |
|||
End |
|||
Do i = 1 To n |
|||
x = cos(pi*(i-0.25)/(n+0.5),prec) |
|||
Do iter = 1 To 10 |
|||
f = p1.1; df = 0 |
|||
Do k = 2 To p1.0 |
|||
df = f + x*df |
|||
f = p1.k + x * f |
|||
End |
|||
dx = f / df |
|||
x = x - dx |
|||
If abs(dx) < epsilon then leave |
|||
End |
|||
r.1.i = x |
|||
r.2.i = 2/((1-x**2)*df**2) |
|||
End |
|||
Return |
|||
cos: Procedure |
|||
/* REXX **************************************************************** |
|||
* Return cos(x) -- with specified precision |
|||
* cos(x) = 1-(x**2/2!)+(x**4/4!)-(x**6/6!)+-... |
|||
* 920903 Walter Pachl |
|||
***********************************************************************/ |
|||
Parse Arg x,prec |
|||
If prec='' Then prec=9 |
|||
Numeric Digits (2*prec) |
|||
Numeric Fuzz 3 |
|||
o=1 |
|||
u=1 |
|||
r=1 |
|||
Do i=1 By 2 |
|||
ra=r |
|||
o=-o*x*x |
|||
u=u*i*(i+1) |
|||
r=r+(o/u) |
|||
If r=ra Then Leave |
|||
End |
|||
Numeric Digits prec |
|||
Return r+0 |
|||
exp: Procedure |
|||
/*********************************************************************** |
|||
* Return exp(x) -- with reasonable precision |
|||
* 920903 Walter Pachl |
|||
***********************************************************************/ |
|||
Parse Arg x,prec |
|||
If prec<9 Then prec=9 |
|||
Numeric Digits (2*prec) |
|||
Numeric Fuzz 3 |
|||
o=1 |
|||
u=1 |
|||
r=1 |
|||
Do i=1 By 1 |
|||
ra=r |
|||
o=o*x |
|||
u=u*i |
|||
r=r+(o/u) |
|||
If r=ra Then Leave |
|||
End |
|||
Numeric Digits (prec) |
|||
Return r+0</lang> |
|||
Output: |
|||
<pre> 1 6.0000000000000000000000000000000000000000 -1.4036E+1 |
|||
2 17.4874646410555689643606840462449458421154 -2.5483 |
|||
3 19.8536919968055821921309108927158495960775 -1.8206E-1 |
|||
4 20.0286883952907008527738054439857661647073 -7.0615E-3 |
|||
5 20.0355777183855621539285357252750939315016 -1.7214E-4 |
|||
6 20.0357469750923438830654575585499253741530 -2.8797E-6 |
|||
7 20.0357498197266007755718729372891903369401 -3.5093E-8 |
|||
8 20.0357498544945172882260918041683132616237 -3.2529E-10 |
|||
9 20.0357498548174338368864419454858704839263 -2.3700E-12 |
|||
10 20.0357498548197898711175766908543458234008 -1.3927E-14 |
|||
11 20.0357498548198037305529147159697031241994 -6.7396E-17 |
|||
12 20.0357498548198037976759531014454017742327 -2.7323E-19 |
|||
13 20.0357498548198037979482458119092690701863 -9.4143E-22 |
|||
14 20.0357498548198037979491844483599375945130 -2.7906E-24 |
|||
15 20.0357498548198037979491872317401917248453 -7.1915E-27 |
|||
16 20.0357498548198037979491872389153958789316 -1.6260E-29 |
|||
17 20.0357498548198037979491872389316236038179 -3.2517E-32 |
|||
18 20.0357498548198037979491872389316560624361 -5.7920E-35 |
|||
19 20.0357498548198037979491872389316561202637 -9.2480E-38 |
|||
20 20.0357498548198037979491872389316561203561 -1.3311E-40 |
|||
20.0357498548198037979491872389316561203562082463657269288113 (exact)</pre> |
|||
=={{header|Tcl}}== |
=={{header|Tcl}}== |