Numerical integration/Gauss-Legendre Quadrature: Difference between revisions
Numerical integration/Gauss-Legendre Quadrature (view source)
Revision as of 14:19, 31 October 2013
, 10 years agoREXX added (translated from PL/I)
Walterpachl (talk | contribs) (→{{header|PL/I}}: correct one declaration for the better!?!) |
Walterpachl (talk | contribs) (REXX added (translated from PL/I)) |
||
Line 989:
</lang>
[[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}}==
| |||